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.

As shown on a previous post, if we define $\theta_1$ and $\theta_2$ as the angles which the first and second rods make with the vertical direction respectively, the following expression can be obtained for the Lagrangian of the system: $$ \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} $$ From the Lagrangian, one can obtain the canonical momenta of the system: $$ \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) \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_ptheta1} \\[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) \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_ptheta2} \end{eqnarray} $$ The Hamiltonian of the system is given by: $$ H = \sum_{i=1}^2 \dot{\theta}_i p_{\theta_i} - L \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian} $$ From the Hamiltonian $H$, we can obtain a set of equations of motion for the system which are equivalent to the Euler-Lagrange equations: $$ \displaystyle \dot{\theta}_i = \frac{\partial H}{\partial p_{\theta_i}}, \quad\quad \dot{p}_{\theta_i} = -\frac{\partial H}{\partial \theta_i} \quad \textrm{ for } \quad i = 1,2 \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eqs} $$ Equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eqs} requires us to write $H$ as a function of the variables $\theta_1$, $\theta_2$, $p_{\theta_1}$ and $p_{\theta_2}$, so from equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian} we see that we must determine $\dot{\theta}_i$ and $L$ in terms of these variables. With that in mind, notice that equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_ptheta1} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_ptheta2} can be written in matrix form as shown below: $$ \left(\begin{matrix} p_{\theta_1} \\[4pt] p_{\theta_1} \end{matrix}\right) = B \left(\begin{matrix} \dot{\theta}_1 \\ \dot{\theta}_2 \end{matrix}\right) \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_p_func_dottheta} $$ where $B$ is a $2 \times 2$ matrix whose entries depend on $\theta_1$ and $\theta_2$: $$ B = \left(\begin{matrix} (m_1 + m_2) l_1^2 & m_2 l_1 l_2 \cos(\theta_1-\theta_2) \\ m_2 l_1 l_2 \cos(\theta_1-\theta_2) & m_2 l_2^2 \end{matrix}\right) $$ From equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_p_func_dottheta}, we can obtain the generalized velocities $\dot{\theta}_i$ in terms of of the canonical momenta $p_{\theta_i}$ and the angles $\theta_i$: $$ \left(\begin{matrix} \dot{\theta}_1 \\ \dot{\theta}_2 \end{matrix}\right) = B^{-1} \left(\begin{matrix} p_{\theta_1} \\[4pt] p_{\theta_1} \end{matrix}\right) \label{%INDEX_eq_dotheta_func_p} $$ The matrix $B$ is indeed invertible for all values of $\theta_1$ and $\theta_2$ since: $$ \begin{eqnarray} \det(B) &=& m_1m_2l_1^2 l_2^2 + m_2^2l_1^2 l_2^2\left[1 - \cos^2(\theta_1 - \theta_2)\right] \nonumber\\[5pt] &=& m_1m_2l_1^2 l_2^2 + m_2^2l_1^2 l_2^2\sin^2(\theta_1 - \theta_2) \nonumber\\[5pt] &\geq& m_1m_2l_1^2 l_2^2 \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_det_B} \end{eqnarray} $$ Being a $2 \times 2$ matrix, $B$ can be inverted directly: $$ B^{-1} = \displaystyle\frac{1}{\det(B)} \left( \begin{matrix} m_2 l_2^2 & -m_2 l_1 l_2 \cos(\theta_1-\theta_2) \\ -m_2 l_1 l_2 \cos(\theta_1-\theta_2) & (m_1 + m_2) l_1^2 \end{matrix} \right) \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_inv_B} $$ Using equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_inv_B} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_det_B} on equation \eqref{%INDEX_eq_dotheta_func_p}, we get (after canceling out common factors and rearranging some terms): $$ \begin{eqnarray} \displaystyle \dot{\theta}_1 &=& \frac{l_2 p_{\theta_1} - l_1 p_{\theta_2}\cos(\theta_1-\theta_2)} {l_1^2l_2\left[m_1 + m_2 \sin^2(\theta_1 - \theta_2)\right]} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta1_func_p}\\[5pt] \displaystyle \dot{\theta}_2 &=& \frac{-m_2 l_2 p_{\theta_1}\cos(\theta_1-\theta_2) + (m_1 + m_2)l_1 p_{\theta_2}} {m_2l_1 l_2^2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta2_func_p} \end{eqnarray} $$ Using equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta1_func_p} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta2_func_p} on equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian} yields the Hamiltonian $H$ in terms of $\theta_1$, $\theta_2$, $p_{\theta_1}$ and $p_{\theta_2}$ (after quite a bit of algebraic work): $$ \begin{eqnarray} H &=& \frac{m_2 l_2^2 p^2_{\theta_1} + (m_1 + m_2)l_1^2 p^2_{\theta_2} - 2m_2 l_1 l_2 p_{\theta_1}p_{\theta_2}\cos(\theta_1-\theta_2)} {2m_2 l_1^2l_2^2 \left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]}\nonumber\\[3pt] & & \quad\quad- (m_1 + m_2) g l_1 \cos\theta_1 - m_2 g l_2\cos\theta_2 \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_right_vars} \end{eqnarray} $$ Equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_right_vars} can finally be used on equation \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eqs} to give us the Hamiltonian equations of motion for the double pendulum: $$ \begin{eqnarray} \displaystyle \dot{\theta}_1 &=& \;\;\frac{\partial H}{\partial p_{\theta_1}} &=& \frac{l_2 p_{\theta_1} - l_1 p_{\theta_2}\cos(\theta_1-\theta_2)} {l_1^2l_2\left[m_1 + m_2 \sin^2(\theta_1 - \theta_2)\right]} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_1} \\[5pt] \displaystyle \dot{\theta}_2 &=& \;\;\frac{\partial H}{\partial p_{\theta_2}} &=& \frac{-m_2 l_2 p_{\theta_1}\cos(\theta_1-\theta_2) + (m_1 + m_2)l_1 p_{\theta_2}} {m_2l_1 l_2^2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_2} \\[5pt] \dot{p}_{\theta_1} &=& -\frac{\partial H}{\partial \theta_1} &=& -(m_1 + m_2)gl_1 \sin\theta_1 - h_1 + h_2\sin\left[2(\theta_1 - \theta_2)\right] \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_3} \\[5pt] \dot{p}_{\theta_2} &=& -\frac{\partial H}{\partial \theta_2} &=& -m_2gl_2 \sin\theta_2 + h_1 - h_2\sin\left[2(\theta_1 - \theta_2)\right] \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_4} \end{eqnarray} $$ where: $$ h_1 = \frac{p_{\theta_1}p_{\theta_2}\sin(\theta_1 - \theta_2)} {l_1l_2\left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]} \\[5pt] \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_h1} $$ $$ h_2 = \frac{m_2 l_2^2 p^2_{\theta_1} + (m_1 + m_2)l_1^2 p^2_{\theta_2} - 2 m_2 l_1l_2 p_{\theta_1}p_{\theta_2}\cos(\theta_1 - \theta_2)} {2l_1^2l_2^2 \left[m_1 + m_2\sin^2(\theta_1 - \theta_2)\right]^2} \label{post_e5ac36fcb129ce95a61f8e8ce0572dbf_h2} $$ Notice that equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_1} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_2} are the same as equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta1_func_p} and \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_eq_dottheta2_func_p}, so in the process of determining the Hamiltonian in terms of the canonical momenta $p_{\theta_i}$ and the angles $\theta_i$, we ended up already obtaining two of the Hamiltonian equations.

Equations \eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_1}-\eqref{post_e5ac36fcb129ce95a61f8e8ce0572dbf_hamiltonian_eq_4} form a set of coupled first order differential equations on the variables $\theta_1$, $\theta_2$, $p_{\theta_1}$ and $p_{\theta_2}$. This set of equations can be solved numerically using a Runge-Kutta (RK) method. A simulator based on the fourth-order RK method can be found here.

References