I also check the GSL (GNU Scientific Library). They have a nice
numerical integrator, but it doesn't allow for a mass matrix.
~Luke
On Sep 29, 4:44 pm, Luke <hazelnu...@gmail.com> wrote:
> Yes, this is something I should look into. I am pretty sure that the
> Netlib codes have this functionality, but it hasn't been wrapped into
> the Python scientific packages that I know of, at least not yet.
> scipy.integrate has odeint and ode, but both need everything in first
> order form, no mass matrix is allowed.
>
> I figured that if I was going to work symbolically, I may as well go
> as far as I can go, so if I can invert the matrices symbolically, I'd
> prefer that. This would also be nice because then the equations could
> be integrated with any numerical integrator out there, as the ode
> formulation would be as generic as possible.
>
> ~Luke
>
> On Sep 29, 4:15 pm, Alan Bromborsky <abro...@verizon.net> wrote:
>
>
>
> > Luke wrote:
>
> > > On Sep 29, 1:09 pm, Ondrej Certik <ond...@certik.cz> wrote:
>
> > >> On Tue, Sep 29, 2009 at 12:49 PM, Luke <hazelnu...@gmail.com> wrote:
>
> > >>> I'm using Sympy from within PyDy to generate the equations of motion for
> > >>> mechanical systems. At the end of the day, the equations can be most
> > >>> generally written as:
> > >>> M(x) * x'' = F(x, x', t)
> > >>> M(x) is what is known as the mass matrix, and will in general depend on
> > >>> the
> > >>> configuration of the system (positions and angles). This matrix needs
> > >>> to be
> > >>> inverted in order to solve for x'', which then will allow for numerical
> > >>> integration or stability analysis.
> > >>> I am generating M(x) symbolically, and in some case is it sparse, but in
> > >>> many cases, it isn't. Each entry of the matrix is a Sympy expression,
> > >>> but
> > >>> instead of trying to invert a matrix of Sympy expressions, I introduce a
> > >>> dummy symbol for each non-zero entries and then invert the matrix of
> > >>> dummy
> > >>> symbols. Humble systems might have 5 degrees of freedom (so a 5x5
> > >>> needs to
> > >>> be inverted), so this inversion isn't so bad, but beyond this, the
> > >>> matrix
> > >>> inversions take a really long time, especially when the matrices are
> > >>> full
> > >>> (no zero entries).
> > >>> I was thinking that it might be nice to have pre-computed matrix
> > >>> inverses
> > >>> for n x n matrices. Matrix inversion is O(n^3), so it would be nice to
> > >>> have
> > >>> all this precomputed symbolically, and this would greatly speed up
> > >>> Sympy's
> > >>> matrix capabilities. Inverses up to say 100x100 could be computed (or
> > >>> maybe
> > >>> something smaller), and then when you need an inverse, everything would
> > >>> be
> > >>> fast. This could also be used behind the scenes (by introduction of
> > >>> symbolic substitution dictionaries) for inverting a matrix full of sympy
> > >>> expressions.
>
> > >> The inversion of a 100x100 dense matrix would be quite a big mess,
> > >> wouldn't it?
>
> > > Yes, it would be disaster.
>
> > >> But I see what you mean, and I think it makes sense to cache it, if it
> > >> speeds things up. But see below.
>
> > >>> Does anybody know if this has been done by somebody somewhere, or have
> > >>> any
> > >>> other ideas on how it could be done better than the way I suggested?
>
> > >> I would first try to profile the inversion code to see why it is slow.
> > >> Because for example the adjoint method is just a ratio of two
> > >> determinants, so it may be the determinants calculation that is slow
> > >> (and should be cached). However, there are 100^2 different
> > >> determinants (right?), but I guess caching 10000 expressions should
> > >> still be doable. But if this is the case, we should just speed up the
> > >> determinant calculation, imho. Looking at the code, we have 2
> > >> algorithms implemented, bareis and berkowitz.
>
> > > I did some timings of the three matrix inversion (GE, ADJ, LU) using
> > > the timeit module. I also timed the two determinant methods as well
> > > to see they the two perform side by side. It seems that bareis (the
> > > default one) is drastically slower than berkowitz, even when you
> > > call .expand() on the berkowitz determinant to make it end up with
> > > identical symbolic expressions. Maybe this should be the default
> > > method for .det()? I went up to 5x5 matrices, it started to get
> > > really slow after that. Gaussian elimination seems to be the slowest
> > > one, at least for dense matrices like the ones I used. Here is the
> > > code I used to do the timings:
>
> > > import timeit
> > > from numpy import zeros, max
> > > import matplotlib.pyplot as plt
>
> > > # Dimension of matrix to invert
> > > n = range(2, 6)
> > > # Number of times to invert
> > > number = 20
> > > # Store the results
> > > t = zeros((len(n), 5))
>
> > > for i in n:
> > > setup_code = "from sympy import Matrix, Symbol\nM = Matrix(%d,\
> > > %d"%(i,i)+",lambda i,j: Symbol('m%d%d'%(i,j)))"
> > > t[i-2, 0] = timeit.Timer('M.inverse_GE()', setup_code).timeit
> > > (number)
> > > t[i-2, 1] = timeit.Timer('M.inverse_ADJ()', setup_code).timeit
> > > (number)
> > > t[i-2, 2] = timeit.Timer('M.inverse_LU()', setup_code).timeit
> > > (number)
> > > t[i-2, 3] = timeit.Timer('M.det(method="bareis")',
> > > setup_code).timeit(number)
> > > t[i-2, 4] = timeit.Timer('M.det(method="berkowitz").expand()',
> > > setup_code).timeit(number)
>
> > > plt.plot(n, t[:,0]/number, label='GE')
> > > plt.plot(n, t[:,1]/number, label='ADJ')
> > > plt.plot(n, t[:,2]/number, label='LU')
> > > plt.plot(n, t[:,3]/number, label='bareis_det')
> > > plt.plot(n, t[:,4]/number, label='berkowitz_det.expand()')
> > > plt.legend(loc=0)
> > > plt.title('Average time to complete 1 matrix inversion/determinant')
> > > plt.xlabel('matrix dimension')
> > > plt.ylabel('Time [seconds]')
> > > plt.xticks(n)
> > > plt.axis([2, n[-1], 0, max(t/number)])
> > > plt.show()
>
> > > I'd be curious to know if others get similar results as I do. I
> > > posted the results of the above script here:
> > >http://i35.tinypic.com/so09hw.jpg
>
> > > It looks like Gaussian elimination is suffering from the bottleneck in
> > > the Bareis determinant since it has an assertion that calls .det() to
> > > make sure the determinant is non-zero.
>
> > >> Also how about the creation of so many expressions, it could be slow too.
>
> > > Yeah, this is true. It seems to me though that if you had all the
> > > inverses computed for fully populated matrices, then if some of the
> > > entries were zero these things would simplify some (or a lot), and
> > > then those expressions would be used rather than the full ones.
>
> > > Maybe it could be an optional module that gets imported in when you
> > > need to work with large matrices and want to have inverses precomputed
> > > so that you don't have to sit and wait for them for hours? That way
> > > it wouldn't affect the normal behavior of Sympy, but could be imported
> > > if desired.
>
> > >> So what I want to say is that we should profile it and determine
> > >> exactly what things are slowing things down, then we should speed them
> > >> up (or cache them).
>
> > > Should I use cProfile for this? Or do you recommend another profiler?
>
> > > ~Luke
>
> > >> Ondrej
>
> > Are there differential equation solvers where you don't have to invert
> > the matrix?
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy" group.
To post to this group, send email to sympy@googlegroups.com
To unsubscribe from this group, send email to sympy+unsubscr...@googlegroups.com
For more options, visit this group at http://groups.google.com/group/sympy?hl=en
-~----------~----~----~----~------~----~------~--~---
