# Posts on Physics

## A proof of Earnshaw's theorem

Posted by Diego Assencio on 2014.09.19 under Physics (Electromagnetism)

Consider a set of $n$ fixed point charges $q_1, q_2, \ldots, q_n$ in vacuum and let ${\bf E}({\bf x})$ be the electric field produced by these charges (see figure 1). Can we find a point ${\bf z}$ in space such that if a charge $q$ is placed there, it will remain in stable equilibrium?

 Fig. 1: A set of $n$ fixed charges and a very small Gaussian surface $S$ centered at a point ${\bf z}$. For ${\bf z}$ to be a point of stable equilibrium (assuming we would place a positive charge $q$ there), ${\bf E}({\bf z})$ must be zero and ${\bf E}({\bf x})$ must point towards ${\bf z}$ at all points ${\bf x}$ of $S$.

The answer is no, and this fact is referred to as the Earnshaw's theorem. We will prove this assuming $q \gt 0$, but the proof is similar for $q \lt 0$.

Let's first think about what do we mean by stable equilibrium for a charge $q$ placed at ${\bf z}$:

 1 the electric force $q{\bf E}({\bf z})$ acting on $q$ must be zero, i.e, ${\bf E}({\bf z}) = {\bf 0}$ 2 for a small displacement $\delta{\bf z}$ of $q$ around ${\bf z}$, $q{\bf E}({\bf z} + \delta{\bf z})$ must point towards ${\bf z}$.

The second property means the electric force on $q$ must be restoring, meaning if we move $q$ a little bit away from ${\bf z}$, this force will will tend to bring $q$ back to ${\bf z}$. Since we are assuming $q \gt 0$, ${\bf E}({\bf z} + \delta{\bf z})$ itself must then point towards ${\bf z}$ for the electric force on $q$ to be restoring.

Consider now a very small Gaussian surface $S$ shaped like a sphere of radius $r$ and centered at ${\bf z}$ (as shown in figure 1). Since ${\bf E}({\bf x})$ points towards $q$ (i.e., inwards) for any point ${\bf x}$ on the surface of this sphere, we have that: $$\int_{S} {\bf E}\cdot d{\bf A} \lt 0 \label{post_bc04395b103021d338b4e30a061bfc74_surface_int}$$ where $d{\bf A}$ is an outward-oriented infinitesimal surface area element. But using the divergence theorem and Gauss's law on equation \eqref{post_bc04395b103021d338b4e30a061bfc74_surface_int} implies we must have a negative charge $Q_{\textrm{enc}}$ enclosed inside our sphere: $$\int_{S} {\bf E}\cdot d{\bf A} = \int_{V}\nabla\cdot{\bf E}\,dV = \displaystyle\frac{Q_{\textrm{enc}}}{\epsilon_0} \lt 0$$ This goes against our initial assumption that ${\bf z}$ is a point in vacuum where no charge is placed ($q$ itself does not count since it does not contribute to the electric field ${\bf E}$).

One of the consequences of Earnshaw's theorem is the fact that a classical model for the atom cannot be stable if it is based on a static (fixed) distribution of charges. For such an atom, the charge configuration would not be in stable equilibrium and therefore unable to resist even the slightest disturbance.

## The double pendulum: Hamiltonian formulation

Posted by Diego Assencio on 2014.05.17 under Physics (Mechanics)

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.

 Fig. 1: A double pendulum.

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.

## The double pendulum: Lagrangian formulation

Posted by Diego Assencio on 2014.02.28 under Physics (Mechanics)

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.

 Fig. 1: A double pendulum.

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_1 \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.