ODE solver from Lagrangian/Variational Methods in C++

Question

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?

Answer 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.