# Convert an implicit ODE system to an explicit ODE set to use Runge-Kutta

For my series on mechanical systems with Lagrange, I would like to add the double pendulum and a rolling pendulum. I set up the Lagrangian $L$ and solved for the ODE of motions, and got the following: $$\ddot x = - \frac l{2m} \cos(\phi) \ddot \phi + \frac l{2m} \sin(\phi) \dot \phi^2$$ $$\ddot \phi = - \frac m{2l} \sin(\phi) \dot \phi + \frac{mg}l \sin(\phi) - \frac{m}{2l} \ddot x \cos(\phi) + \frac{m}{2l} \dot x \sin(\phi) \dot \phi$$

I am using the ODE solver scipy.integrate.odeint, and it takes a function $\vec f$ so that: $$\frac{\mathrm d}{\mathrm dt} \vec y(t) = \vec f\left((\vec y(t), t\right)$$

If the functions were not implicit, I would just do the following approach: $$\frac{\mathrm d}{\mathrm dt} \begin{pmatrix} x \\ \phi \\ \dot x \\ \dot \phi \end{pmatrix} = \begin{pmatrix} \dot x \\ \dot \phi \\ \text{RHS of first equation} \\ \text{RHS of second equation} \end{pmatrix}$$

My problem is that I do not know $\ddot x$ and $\ddot\phi$ to plug in into the the calculations. So far, I have tried putting the second (hand written) equation into the first, then calculating $\ddot x$. Then I would calculate $\ddot \phi$ with that value of $\ddot x$. The result were strange results and stability issues.

How would I integrate this entangled system of equations?

# Update 2013-09-22 14:56:21+02:00

My current function $\vec f$ is the following. I plugged $\ddot \phi$, i. e. the second equation, into the first equation and solved for $\ddot x$ so that I have an equation in $\ddot x$ without any other second derivatives on the RHS. My Python code is the following, where I just divided by $1 - \cos^2(\phi)/4$ in the first equation:

x = y[0]
phi = y[1]
dot_x = y[2]
dot_phi = y[3]

ddt_x = dot_x
ddt_phi = dot_phi

ddt_dot_x = (- self.l / (2 * self.m) * np.cos(phi) * (- self.m / (2 * self.l) * np.sin(phi) * dot_phi + self.m * self.g / self.l * np.sin(phi) + self.m / (2 * self.l) * dot_x * np.sin(phi) * dot_phi) + self.l / (2 * self.m) * np.sin(phi) * dot_phi**2) / (1 - 1/4 * np.cos(phi)**2)

ddt_dot_phi = - self.m / (2 * self.l) * np.sin(phi) * dot_phi + self.m * self.g / self.l * np.sin(phi) - self.m / (2 * self.l) * ddt_dot_x * np.sin(phi) + self.m / (2 * self.l) * dot_x * np.sin(phi) * dot_phi

return [ddt_x, ddt_phi, ddt_dot_x, ddt_dot_phi]


When I try to integrate this, I get results that do not really reflect any physical effect. You can see that in the animation that does a sudden twitch.

• Are you certain that you have $\ddot\phi$ appearing also on the right hand side? I'm not saying that you don't, but it's definitely a case that one doesn't usually see. – Wolfgang Bangerth Sep 21 '13 at 2:40
• @WolfgangBangerth I think so. You can check my derivation if you can spare the time. It is just the first page of the PDF. – Martin Ueding Sep 21 '13 at 10:49
• It doesn't look obviously wrong to me. I think you should make the substitution at the top of page 2, so that you get an equation where you have $(1-[\cos \phi]^2/4)\ddot x$ on the left and only first derivatives on the right. Note that the pre-factor is nicely behaved. Then use the original second equation, which has the same structure. You can easily convert this into a first order system of standard form. – Wolfgang Bangerth Sep 21 '13 at 11:09
• @WolfgangBangerth I tried that, see my edits in the question, but that does not seem to give a good solution. – Martin Ueding Sep 22 '13 at 13:04
• You'll have to debug your code, I'm afraid. – Wolfgang Bangerth Sep 22 '13 at 21:31

My first thought would be to convert this system to a first-order system of differential algebraic equations.

In essence, you'd set something like:

\begin{align} y &= \dot{x}, \\ \dot{y} &= \ddot{x}, \\ \gamma &= \dot{\varphi}, \\ \dot{\gamma} &= \ddot{\varphi} \end{align}

and then rearrange your implicit ODE so it looks something like:

\begin{align} f(x, \dot{x}, y, \dot{y}, \varphi, \dot{\varphi}, \gamma, \dot{\gamma}) = 0, \end{align}

which is a DAE.

The potential pitfall of this approach is that your resulting DAE system might be high-index. You should see if the DAE system can be expressed as a Hessenberg index-1 or Hessenberg index-2 DAE, for which there good numerical methods available (for instance, IDA in the SUNDIALS package, accessed in Python from scikits.odes, from Assimulo, or another Python interface to IDA). If it can't, there will likely be stability issues in calculating a solution. Most methods are geared towards solving problems with an index of less than 3. In particular, mechanical systems are known for high-index formulations. (See the book by Ascher and Petzold for some good examples.) There are many definitions of index, so I'm just going to point you to the Scholarpedia article on DAEs for details.

There are also DAE solvers in PETSc accessible through petsc4py, but it's not clear to me from the documentation if they can solve problems with index greater than 1.

• If you look at the wikipedia article about the double pendulum, you see that the four last equations are matching my approach. I am wondering how they got that done. The two equations that I gave are for a rolling pendulum, but the problem is the same, just a little simpler. – Martin Ueding Sep 20 '13 at 18:18
• My guess is that they took the two equations defining the angular momenta and solved that system to express the time derivatives of the angles (angular velocities?) in terms of the momenta. Then, it seems like the differential equations for the angular momenta follow from the Lagrangian. (I've only taken a class on introductory classical mechanics, so that step isn't so clear to me.) That formulation should be preferable for numerical integration, but without its spectral properties or conservation relationships, I couldn't recommend a specific method. – Geoff Oxberry Sep 20 '13 at 20:26
• Probably the biggest challenge integrating higher-index DAE is producing consistent initial conditions. If you have consistent initial conditions, many of the PETSc solvers will work out of the box, but we don't have special support for higher index initial conditions. – Jed Brown Sep 22 '13 at 13:53

You certainly have a number of cursoriness mistakes in your derivation, but often only part of the mistake is propagated to later equations.

1. $V=mgly_2$ should be $V=mgy_2$.
2. $L=\frac{m}{2}\left(2\dot x^2+\dot x l\cos(\phi)\dot\phi\right)-mgl\cos(\phi)$ should be $L=\frac{m}{2}\left(2\dot x^2+2\dot x l\cos(\phi)\dot\phi\right)-mgl\cos(\phi)$.
3. $\pi_x=\frac{\partial L}{\partial\dot x}=2m\dot x+l\cos(\phi)\dot\phi$ should be $\pi_x=\frac{\partial L}{\partial\dot x}=2m\dot x+ml\cos(\phi)\dot\phi$

The following equations are already tainted by these mistakes, so this should be enough for the moment.

Regarding DAEs and their index, the Lagrangian formulation for your problem is uncritical, and leads to a system of index $\leq 1$. So I can assure you that the main problem are the cursoriness mistakes all over the place.