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:00My 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_xddt_phi = dot_phiddt_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_phireturn [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.
Convert an implicit ODE system to an explicit ODE set to use Runge-Kutta
ode;runge kutta