5

I have a general question which I will phrase in the context of a more concrete situation.

If one wants to find the dynamics of a double pendulum, one can mathematically derive the equations of motions, rewrite the ODEs to have a special form useful for numerical computation and solve the ODEs using, say, odeint in c++ (see an example of this on stack overflow https://stackoverflow.com/a/30582741).

Now imagine we want to do the same for n coupled pendula (n known at run time). This requires us to write a so-called Lagrangian (kinetic energy - potential energy) and different derivatives of this functions will be the ODEs we need to solve. Furthermore these ODEs must be rewritten in a form suitable for odeint. This can for general n be hard to do by hand.

In programs like Mathematica and Maple, this is actually quite easy. One can explicitly derive the ODEs needed from the Lagrangian and the ODE solver do not need us to put the equations in any special form (see example in mathematica here https://mathematica.stackexchange.com/a/84279).

Is it possible to do such a thing in c++ without going through too much trouble?

Potential Approach:

A potential approach could be to use the c++ package ginac. This can help us to analytically derive the ODEs. But I do not know how to rewrite the expressions coming from ginac, into a form suitable for numerical computation in odeint. Any ideas?

Community
  • 1
  • 1
Heidar
  • 159
  • 1
  • 1
    For n known only at run time, I'm afraid you are facing almost impossible difficulties. If you know n at compile time I think you could achieve such a computation based on automatic differentiation and expression templates. From my (superficial) check of ginac I think this library is also based on expression templates, hence would need n at compile time. – mariomulansky Jul 07 '15 at 11:57
  • @mariomulansky Thanks a lot for the comment. For now I would be happy if I can do it in compile time. How can I use expression templates to convert the ginac output into a form odeint understands? Does there exist some simple standard way? – Heidar Jul 08 '15 at 00:22
  • 1
    I have no experience with GiNaC, so I can't be of much help. However, reading more through their tutorial, it occurs to me that you could be able to construct an expression at runtime, differentiate it and evaluate the result. So my first comment above might be incorrect. Trying to write the Langragian as a GiNaC expression, differentiate it and use the resulting expressions in the rhs for odeint sounds does sound like a feasible approach. If you do get this running, I would be very interested in the code. Maybe you can share it on odeint's github: https://github.com/headmyshoulder/odeint-v2 – mariomulansky Jul 08 '15 at 08:59
  • 1
    You could also treat the Lagrange system as an index-3 DAE system and use the BDF or RADAU solvers that are commonly used for constraint dynamics. Or manually differentiate the constraint equations twice and numerically solve the resulting non-linear system for the second derivatives via Newton in the ODE function. – Lutz Lehmann Jul 08 '15 at 16:53
  • @mariomulansky Yeah it seems that one can directly evaluate the expressions, which should be enough to get it into odeint. I will try to implement it when I get some time, I'll be happy to share the code. For this more specific case (n-pendula) there might be a simpler approach. If I use normal cartesian coordinates and implement the pendula constraints using lagrange multipliers, the Hamiltonian will be of the form H=T(p) + V(q) (multipliers among q's). Maybe from there one can directly use odeint's symplectic integrators? – Heidar Jul 09 '15 at 01:21
  • @LutzL Thanks for the suggestions. I will try to read a bit about DAE systems and the BDF/RADAU solvers. Could be interesting to implement it like that. – Heidar Jul 09 '15 at 01:24

1 Answers1

1

The trivial intertia term helps in the sense that you only need to compute dV/dq and not dT/dp. odeint provides a version of the symplectic integrator that only expects dV/dq and assumes dT/dp being trivial as in your case. However, you still need to obtain a derivative.

mariomulansky
  • 963
  • 5
  • 7