In the formulation I use for PyDy, the equations of motion are
generated in first order form. For holonomic systems with n degrees
of freedom, there will be 2n first order equations and the first n of
these I refer to as the kinematic differential equations. In the
simplest case, the form of the first n equations would be like:
q_i' = u_i (i = 1,...,n)
The remaining n ODE's (dynamic differential equations) would be the
ones with the mass matrix in the form like:
M(q)*u' = F(q, u, t)
Where q, u are n x 1 vectors, M is nxn, and F: n x n x 1 --> n
Usually this choice of 'generalized speeds' (a.k.a quasi-coordinates)
which result in kinematic differential equations like: "q_i' = u_i" is
not the best choice out there. But even with other choices, they can
be written as:
q' = T(q) * u + R(t)
Where in most cases R(t) is zero, except when you have some prescribed
motion in your system. The dynamic differential equations would still
have the same form as before.
So they are linear in the generalized speeds, but the T(q) matrix
would in general be comprised of trigonometric entries which depend on
the configuration and/or parameters. Such choices of generalized
speeds can often lead to reduced complexity of the dynamic equations,
with some (usually worthwhile) trade-off in the complexity of the
kinematic equations.
For nonholonomic systems, depending on how you formulate things you
may end up having to integrate fewer dynamic equations. For example,
in the rolling disc, there are 5 coordinates that configure the disc,
but only 3 degrees of freedom. If you are interested in all of the
motion, you would integrate 5 kinematic D.E.'s and 3 dynamic D.E.'s.
If you only are interested in the dynamic equations, you would end up
only integrating 1 kinematic D.E. (the one associated with the lean of
the disc), and 3 dynamic D.E.'s. So in these types of systems, it
depends on what your needs are.
~Luke
On Sep 29, 4:22 pm, Alan Bromborsky <abro...@verizon.net> wrote:
> Alan Bromborsky 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?
>
> Also, what does the system look like if written as a first order system
> of 2n equations?
--~--~---------~--~----~------------~-------~--~----~
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
-~----------~----~----~----~------~----~------~--~---
