Luke wrote:
> 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?
>>>
> >
>
>
Would you give an example set of equations you wish to solve?
--~--~---------~--~----~------------~-------~--~----~
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
-~----------~----~----~----~------~----~------~--~---
