Simplistic simulation of double pendulum doesn't work. Is this approach valid?

Question

Introduction

I'd like to preface this by saying that while I'm a reasonably competent developer, my calculus skills leave a lot to be desired.

I've been reading up on basic physics, and I thought I had a grip on the behaviour of the Lagrangian formulation of physics, but it turns out there is something I just don't understand, which is why I'm turning to this place for assistance.

I've been simulating various systems like pendulums, gravitational systems, etc and my very simplistic approach has worked fine. However, when I tried it on the double pendulum, I'm not getting the expected results.

Calculating the Lagrangian

I'm assuming the coordinates of the masses:

$$ \begin{align} x_1 &= sin(θ_1)\,l_1 \\ x_2 &= x_1 + sin(θ_2)\,l_2 \\ y_1 &= -cos(θ_1)\,l_1 \\ y_2 &= y_1 - cos(θ_2)\,l_2 \end{align} $$

I then use Maxima to compute the Lagrangian, which is:

$$ {{2\,l_{1}\,l_{2}\,m_{2}\,\dot{θ}_{1}\,\dot{θ}_{2}\,\sin θ_{1}\,\sin θ_{2}+ \left(2\,l_{1}\,l_{2}\,m_{2}\,\dot{θ}_{1}\,\dot{θ}_{2}\,\cos θ_{1}+2\,g\,l_{2}\, m_{2}\right)\,\cos θ_{2}+\left(2\,g\,l_{1}\,m_{2}+2\,g\,l_{1}\,m_{1} \right)\,\cos θ_{1}+l_{2}^2\,m_{2}\,\dot{θ}_{2}^2+\left(l_{1}^2\,m_{2}+ l_{1}^2\,m_{1}\right)\,\dot{θ}_{1}^2}\over{2}} $$

I then use the Euler-Lagrange equation to give me the solution for the first coordinate:

$$ \ddot{θ}_1 = -{{\left(\ddot{θ}_{2}\,l_{2}\,m_{2}\,\sin θ_{1}-l_{2}\,m_{2}\,\dot{θ}_{2}^2\, \cos θ_{1}\right)\,\sin θ_{2}+\left(l_{2}\,m_{2}\,\dot{θ}_{2}^2\,\sin θ_{1}+\ddot{θ}_{2}\,l_{2}\,m_{2}\,\cos θ_{1}\right)\,\cos θ_{2}+\left(g\, m_{2}+g\,m_{1}\right)\,\sin θ_{1}}\over{l_{1}\,m_{2}+l_{1}\,m_{1}}} $$

And the solution for the second coordinate:

$$ \ddot{θ}_2 = -{{\left(\ddot{θ}_{1}\,l_{1}\,\sin θ_{1}+l_{1}\,\dot{θ}_{1}^2\,\cos θ_{1}+g \right)\,\sin θ_{2}+\left(\ddot{θ}_{1}\,l_{1}\,\cos θ_{1}-l_{1}\,\dot{θ}_{1}^2\, \sin θ_{1}\right)\,\cos θ_{2}}\over{l_{2}}} $$

While my Lagrangian looks a bit different from the typical one I can find when searching online for the Lagrangian for a double pendulum, they do seem to be equivalent, and I did try to enter the one I found online and I get the same numbers from that one.

Trying to simulate the motion

Once I got this, I wanted to simulate it in the same simplistic manner as I have have simulated other physical systems, by accumulating the positions and velocities.

For each iteration, I compute 6 variables like so:

$$ \texttt{thetadotdot_1 = } \ddot{\theta}_1 \Delta t \texttt{ /* The result above */} \\ \texttt{thetadotdot_2 = } \ddot{\theta}_2 \Delta t \texttt{ /* The result above */} \\ \texttt{thetadot_1 = thetadot_1 + thetadotdot_1} \\ \texttt{thetadot_2 = thetadot_2 + thetadotdot_2} \\ \texttt{theta_1 = theta_1 + thetadot_1} \\ \texttt{theta_2 = theta_2 + thetadot_1} $$

I then draw the two objects based on theta_1 and theta_2 alogn with the values of $l_1$ and $l_2$.

The problem I am facing is that while the two objects swings back and forth, they do so completely independently. Each one swings like a regular pendulum, without being influenced by the other. The fact that they swing at all suggests that the way I approach the problem isn't entirely wrong, but I clearly made a mistake somewhere.

I'd be very happy is someone could point out where my mistake is?

Answer 1

The Lagrangian has the traditional form $L=T-V$ of difference of kinetic and potential energy. $\newcommand{\pd}[2]{\frac{\partial#1}{\partial#2}}$ \begin{align} T&=\frac{m_1}{2}(\dot x_1^2+\dot y_1^2)+\frac{m_2}{2}(\dot x_2^2+\dot y_2^2) \\ &=\frac{m_1}{2}l_1^2\dot\theta_1^2+\frac{m_2}{2}\Bigl(l_1^2\dot\theta_1^2+2l_1l_2\cos(\theta_1-\theta_2)\dot\theta_1\dot\theta_2+l_2^2\dot\theta_2^2\Bigr) \\ V&=m_1gy_1+m_2gy_2 \\ &=-(m_1+m_2)l_1g\cos\theta_1-m_2l_2\cos\theta_2 \end{align} with partial derivatives \begin{align} \pd{L}{\dot\theta_1}=\pd{T}{\dot\theta_1} &=(m_1+m_2)l_1^2\dot\theta_1+m_2l_1l_2\cos(\theta_1-\theta_2)\dot\theta_2 \\ \pd{L}{\dot\theta_2}=\pd{T}{\dot\theta_2} &=m_2l_1l_2\cos(\theta_1-\theta_2)\dot\theta_1+m_2l_2^2\dot\theta_2 \\ \pd{L}{\theta_1}&=-m_2l_1l_2\sin(\theta_1-\theta_2)\dot\theta_1\dot\theta_2-(m_1+m_2)l_1g\sin\theta_1 \\ \pd{L}{\theta_2}&=m_2l_1l_2\sin(\theta_1-\theta_2)\dot\theta_1\dot\theta_2-m_2l_2g\sin\theta_2 \end{align} and Euler-Lagrange equations \begin{align} \frac{d}{dt}\pd{L}{\dot\theta_1}=\pd{L}{\theta_1}&\implies \\ (m_1+m_2)l_1^2\ddot\theta_1+m_2l_1l_2\cos(\theta_1-\theta_2)\ddot\theta_2 &=-m_2l_1l_2\sin(\theta_1-\theta_2)\dot\theta_2^2-(m_1+m_2)l_1g\sin\theta_1 \\ \frac{d}{dt}\pd{L}{\dot\theta_2}=\pd{L}{\theta_2}&\implies \\ m_2l_1l_2\cos(\theta_1-\theta_2)\ddot\theta_1+m_2l_2^2\ddot\theta_2 &=m_2l_1l_2\sin(\theta_1-\theta_2)\dot\theta_1^2-m_2l_2g\sin\theta_2 \end{align}

This now is a nice linear system for the second derivatives $$ \begin{bmatrix} (m_1+m_2)&m_2\cos(\theta_1-\theta_2)\\ m_2\cos(\theta_1-\theta_2)&m_2 \end{bmatrix} \begin{bmatrix} l_1\ddot\theta_1\\l_2\ddot\theta_2 \end{bmatrix} = \begin{bmatrix} -m_2l_2\sin(\theta_1-\theta_2)\dot\theta_2^2-(m_1+m_2)g\sin\theta_1 \\ m_2l_1\sin(\theta_1-\theta_2)\dot\theta_1^2-m_2g\sin\theta_2 \end{bmatrix} $$ or after division by $m_2$, setting $\mu=1+\frac{m_1}{m_2}$, $$ \begin{bmatrix} \mu&\cos(\theta_1-\theta_2)\\ \cos(\theta_1-\theta_2)&1 \end{bmatrix} \begin{bmatrix} l_1\ddot\theta_1\\l_2\ddot\theta_2 \end{bmatrix} =- \begin{bmatrix} \mu g\sin\theta_1+l_2\sin(\theta_1-\theta_2)\dot\theta_2^2 \\ g\sin\theta_2-l_1\sin(\theta_1-\theta_2)\dot\theta_1^2 \end{bmatrix} $$ This can be solved via manual elimination, the Cramer formula/multiplication with the adjoint matrix, or a linear system solver.