Consider the double pendulum shown on figure 1. A double pendulum is formed by attaching a pendulum directly to another one. Each pendulum consists of a bob connected to a massless rigid rod which is only allowed to move along a vertical plane. The pivot of the first pendulum is fixed to a point $O$. All motion is frictionless.

In our discussion, the fixed point $O$ will be taken as the origin of the Cartesian coordinate system with the $x$ axis pointing along the horizontal direction and the $y$ axis pointing vertically upwards. Let $\theta_1$ and $\theta_2$ be the angles which the first and second rods make with the vertical direction respectively. As can be seen on figure 1, the positions of the bobs are given by: $$ \begin{eqnarray} x_1 &=& l_1 \sin\theta_1 \quad\quad & y_1 &=& -l_1 \cos\theta_1\\[5pt] x_2 &=& l_1 \sin\theta_1 + l_2 \sin\theta_2 \quad\quad & y_2 &=& -l_1\cos\theta_1 -l_2\cos\theta_2 \end{eqnarray} $$ Differentiating the quantities above with respect to time, we obtain the velocities of the bobs: $$ \begin{eqnarray} \dot{x}_1 &=& l_1 \dot{\theta}_1\cos\theta_1 \quad\quad & \dot{y}_1 &=& l_1 \dot{\theta}_1\sin\theta_1\\[5pt] \dot{x}_2 &=& l_1 \dot{\theta}_1\cos\theta_1 + l_2 \dot{\theta}_2\cos\theta_2 \quad\quad & \dot{y}_2 &=& l_1 \dot{\theta}_1\sin\theta_1 + l_2\dot{\theta}_2\sin\theta_2 \end{eqnarray} $$

The double pendulum is a very interesting system as it is very simple but can show chaotic behavior for certain initial conditions. In this regime, slightly changing the initial values of the angles ($\theta_1,\theta_2$) and angular velocities ($\dot{\theta}_1,\dot{\theta}_2$) makes the trajectories of the bobs become very different from the original ones.

The Lagrangian for the double pendulum is given by $L = T - V$, where $T$ and $V$ are the kinetic and potential energies of the system respectively. The kinetic energy $T$ is given by: $$ \begin{eqnarray} T &=& \displaystyle\frac{1}{2}m_1v_1^2 + \frac{1}{2}m_2v_2^2 \nonumber \\[5pt] &=& \frac{1}{2}m_1(\dot{x}_1^2 + \dot{y}_1^2) + \frac{1}{2}m_2(\dot{x}_2^2 + \dot{y}_2^2) \nonumber \\[5pt] &=& \frac{1}{2}m_1 l_1^2 \dot{\theta}_1^2 + \frac{1}{2}m_2\left[l_1^2 \dot{\theta}_1^2 + l_2^2 \dot{\theta}_2^2 + 2l_1l_2\dot{\theta}_1\dot{\theta}_2 \cos(\theta_1 - \theta_2)\right] \end{eqnarray} $$ where above we used the fact that $\cos\theta_1\cos\theta_2 + \sin\theta_1\sin\theta_2 =$ $\cos(\theta_1-\theta_2)$. The potential energy $V$ is given by: $$ \begin{eqnarray} V &=& m_1 g y_1 + m_2gy_2 \nonumber \\[5pt] &=& -m_1 g l_1 \cos\theta_1 - m_2 g (l_1 \cos\theta_1 + l_2 \cos\theta_2) \nonumber\\[5pt] &=& -(m_1 + m_2) g l_1 \cos\theta_1 - m_2 g l_2\cos\theta_2 \end{eqnarray} $$ The Lagrangian of the system is then: $$ \begin{eqnarray} L =\frac{1}{2}(m_1 + m_2) l_1^2 \dot{\theta}_1^2 &+& \frac{1}{2}m_2 l_2^2 \dot{\theta}_2^2 + m_2l_1l_2\dot{\theta}_1\dot{\theta}_2 \cos(\theta_1 - \theta_2)\nonumber\\[3pt] &+&(m_1 + m_2) g l_1 \cos\theta_1 + m_2 g l_2\cos\theta_2 \end{eqnarray} $$ The canonical momenta associated with the coordinates $\theta_1$ and $\theta_2$ can be obtained directly from $L$: $$ \begin{eqnarray} \displaystyle p_{\theta_1} &=& \frac{\partial L}{\partial \dot{\theta}_1} = (m_1 + m_2) l_1^2 \dot{\theta}_1 + m_2 l_1 l_2 \dot{\theta}_2\cos(\theta_1-\theta_2) \\[5pt] \displaystyle p_{\theta_2} &=& \frac{\partial L}{\partial \dot{\theta}_2} = m_2 l_2^2 \dot{\theta}_2 + m_2 l_1 l_2 \dot{\theta}_1\cos(\theta_1-\theta_2) \end{eqnarray} $$ The equations of motion of the system are the Euler-Lagrange equations: $$ \displaystyle \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\theta}_i}\right) - \frac{\partial L}{\partial \theta_i} = 0 \Longrightarrow \frac{d p_{\theta_i}}{dt} - \frac{\partial L}{\partial \theta_i} = 0 \quad \textrm{ for } \quad i = 1,2 \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq} $$ Since: $$ \begin{eqnarray} \displaystyle \frac{d p_{\theta_1}}{dt} &=& (m_1 + m_2) l_1^2 \ddot{\theta}_1 + m_2 l_1 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2) \nonumber \\ &\;&\quad -\, m_2 l_1 l_2 \dot{\theta}_2\dot{\theta}_1\sin(\theta_1 - \theta_2) + m_2 l_1 l_2 \dot{\theta}_2^2\sin(\theta_1 - \theta_2) \\[5pt] \displaystyle \frac{d p_{\theta_2}}{dt} &=& m_2 l_2^2 \ddot{\theta}_2 + m_2 l_1 l_2 \ddot{\theta}_1 \cos(\theta_1 - \theta_2) \nonumber \\ &\;&\quad -\, m_2 l_1 l_2 \dot{\theta}_1^2 \sin(\theta_1 - \theta_2) + m_2 l_1 l_2 \dot{\theta}_1\dot{\theta}_2 \sin(\theta_1 - \theta_2)\\[5pt] \displaystyle \frac{\partial L}{\partial \theta_1} &=& - m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1-\theta_2) - (m_1 + m_2) g l_1 \sin\theta_1 \\[5pt] \displaystyle \frac{\partial L}{\partial \theta_2} &=& m_2 l_1 l_2 \dot{\theta}_1 \dot{\theta}_2 \sin(\theta_1-\theta_2) - m_2 g l_2 \sin\theta_2 \end{eqnarray} $$ then equation \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq} yields (after dividing by $l_1$ when $i = 1$ and by $m_2 l_2$ when $i = 2$): $$ \begin{eqnarray} (m_1 + m_2) l_1 \ddot{\theta}_1 &+& m_2 l_2 \ddot{\theta}_2 \cos(\theta_1 - \theta_2) \nonumber\\ &+& m_2 l_2 \dot{\theta}_2^2\sin(\theta_1 - \theta_2) + (m_1 + m_2) g \sin\theta_1 = 0 \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1} \end{eqnarray} $$ $$ l_2 \ddot{\theta}_2 + l_1 \ddot{\theta}_1 \cos(\theta_1 - \theta_2) - l_1 \dot{\theta}_1^2 \sin(\theta_1 - \theta_2) + g \sin\theta_2 = 0 \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2} $$ Equations \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1} and \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2} form a system of coupled second-order nonlinear differential equations. Dividing equation \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1} by $(m_1 + m_2)l_1$ and equation \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2} by $l_2$ and also moving all terms which do not involve $\ddot{\theta}_1$ and $\ddot{\theta}_2$ to the right-hand side, we obtain: $$ \begin{eqnarray} \ddot{\theta}_1 + \alpha_1(\theta_1,\theta_2) \ddot{\theta}_2 &=& f_1(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1_simple}\\[5pt] \ddot{\theta}_2 + \alpha_2(\theta_1,\theta_2) \ddot{\theta}_1 &=& f_2(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) \label{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2_simple} \end{eqnarray} $$ where: $$ \begin{eqnarray} \alpha_1(\theta_1,\theta_2) &:=& \displaystyle\frac{l_2}{l_1}\left(\frac{m_2}{m_1 + m_2}\right)\cos(\theta_1 - \theta_2) \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_alpha1}\\[5pt] \alpha_2(\theta_1,\theta_2) &:=& \frac{l_1}{l_2}\cos(\theta_1-\theta_2) \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_alpha2} \end{eqnarray} $$ and: $$ \begin{eqnarray} \displaystyle f_1(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) &:=& -\frac{l_2}{l_1}\left(\frac{m_2}{m_1+m_2}\right) \dot{\theta}_2^2\sin(\theta_1 - \theta_2) - \frac{g}{l_1} \sin\theta_1 \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_f_1}\\[5pt] \displaystyle f_2(\theta_1, \theta_2, \dot{\theta}_1, \dot{\theta}_2) &:=& \frac{l_1}{l_2}\dot{\theta}_1^2\sin(\theta_1-\theta_2) - \frac{g}{l_2} \sin\theta_2 \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_f_2} \end{eqnarray} $$ Interestingly, $f_1$ does not depent on $\dot{\theta}_1$ and $f_2$ does not depend on $\dot{\theta}_2$. Equations \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq1_simple} and \eqref{post_1500c66ae7ab27bb0106467c68feebc6_euler_lagrange_eq2_simple} can be combined into a single equation: $$ A \left( \begin{matrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{matrix} \right) = \left( \begin{matrix} 1 & \alpha_1 \\[3pt] \alpha_2 & 1 \end{matrix} \right) \left( \begin{matrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{matrix} \right) = \left( \begin{matrix} f_1 \\ f_2 \end{matrix} \right) \label{%INDEX_eq_lin_sys_ddtheta} $$ where the matrix $A$ depends on $\theta_1$ and $\theta_2$ since $\alpha_1$ and $\alpha_2$ depend on these variables. Being a $2 \times 2$ matrix, $A$ can be inverted directly: $$ A^{-1} = \displaystyle\frac{1}{\det(A)} \left( \begin{matrix} 1 & -\alpha_1 \\ -\alpha_2 & 1 \end{matrix} \right) = \displaystyle\frac{1}{1 - \alpha_1\alpha_2} \left( \begin{matrix} 1 & -\alpha_1 \\ -\alpha_2 & 1 \end{matrix} \right) \label{post_1500c66ae7ab27bb0106467c68feebc6_inv_A} $$ Before we continue, notice that $A$ is always invertible since: $$ \det(A) = 1 - \alpha_1\alpha_2 = 1 - \left(\frac{m_2}{m_1 + m_2}\right)\cos^2(\theta_1 - \theta_2) \gt 0 $$ because $m_2 / (m_1 + m_2) \lt 1$ and $\cos^2(x) \leq 1$ for all real values of $x$. From equations \eqref{%INDEX_eq_lin_sys_ddtheta} and \eqref{post_1500c66ae7ab27bb0106467c68feebc6_inv_A} we obtain: $$ \left( \begin{matrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \end{matrix} \right) = A^{-1} \left( \begin{matrix} f_1 \\[3pt] f_2 \end{matrix} \right) = \displaystyle\frac{1}{1 - \alpha_1\alpha_2} \left( \begin{matrix} f_1 - \alpha_1 f_2\\[3pt] -\alpha_2 f_1 + f_2 \end{matrix} \right) $$ Finally, letting $\omega_1 := \dot{\theta}_1$ and $\omega_2 := \dot{\theta}_2$, we can write the equations of motion of the double pendulum as a system of coupled first order differential equations on the variables $\theta_1$, $\theta_2$, $\omega_1$, $\omega_2$: $$ \displaystyle\frac{d}{dt} \left( \begin{matrix} \theta_1 \\[1pt] \theta_2 \\[1pt] \omega_1 \\[1pt] \omega_2 \end{matrix} \right) = \left( \begin{matrix} \omega_1 \\ \omega_2 \\ g_1(\theta_1,\theta_2,\omega_1,\omega_2) \\ g_2(\theta_1,\theta_2,\omega_1,\omega_2) \end{matrix} \right) \label{post_1500c66ae7ab27bb0106467c68feebc6_first_order_eq_theta_omega} $$ $$ g_1 := \displaystyle\frac{f_1 - \alpha_1 f_2}{1 - \alpha_1\alpha_2} \quad\quad g_2 := \displaystyle\frac{-\alpha_2 f_1 + f_2}{1 - \alpha_1\alpha_2} \label{post_1500c66ae7ab27bb0106467c68feebc6_eq_g1_g2} $$ where $\alpha_i = \alpha_i(\theta_1,\theta_2)$ and $f_i = f_i(\theta_1,\theta_2,\omega_1,\omega_2)$ for $i = 1,2$ are given on equations \eqref{post_1500c66ae7ab27bb0106467c68feebc6_eq_alpha1}-\eqref{post_1500c66ae7ab27bb0106467c68feebc6_eq_f_2}.

Equation \eqref{post_1500c66ae7ab27bb0106467c68feebc6_first_order_eq_theta_omega} can be solved numerically using a Runge-Kutta (RK) method. A simulator based on the fourth-order RK method can be found here.

References