Introduction
The double pendulum, or chaotic pendulum, is a system that comprises
of a pendulum attached to the end of another pendulum. This system, as I will derive below,
is governed by a set of coupled non-linear ordinary differential equations. Beyond small angles,
in which case the system can be described as a coupled oscillator, the system cannot be solved analytically
and numerical methods must be employed. The system exhibits properties of non-linear systems
such as chaotic motion and sensitivity to initial conditions.
You can start and stop the simulation with the start/pause button. To adjust the parameters of the system,
the time rate, time step, and path settings, click the settings button to open up the settings menu.
You can also drag the pendulum bobs to modify their angles. Choosing different presets from
the dropdown initializes the system with a set of initial conditions. The reset button resets
the system to whatever initial conditions are selected from the dropdown. Clicking on the graphs tab
gives you options for graphing parameters to create phase-portraits.
Equations of Motion
The double pendulum consists of a mass m₁ suspended by a massless rod of length L₁,
from a fixed pivot, and a second mass m₂ suspended by a massless rod of length L₂ from m₁.
The equations of motion for the system will be derived using the Lagrangian formalism.
We will use θ₁, the angle of m₁ relative to it's vertical,
and θ₂, the angle of m₂ relative to it's vertical, as our generalized coordinates.
The resulting kinetic and potential energy of the system is:
$$T = \frac{1}{2}ML_1^{2}\dot{\theta_1}^{2} + \frac{1}{2}m_2L_2^{2}\dot{\theta_2}^{2} + m_2L_1L_2\dot{\theta_1}\dot{\theta_2}\cos\Delta\theta$$
$$U = -MgL_1\cos\theta_1 - m_2gL_2\cos\theta_2$$
where I have defined M = m₁ + m₂ and Δθ = θ₁ - θ₂.
The resulting Lagrangian is just T - U:
$$\mathcal{L} = \frac{1}{2}ML_1^{2}\dot{\theta_1}^{2} + \frac{1}{2}m_2L_2^{2}\dot{\theta_2}^{2} + m_2L_1L_2\dot{\theta_1}\dot{\theta_2}\cos\Delta\theta + MgL_1\cos\theta_1 + m_2gL_2\cos\theta_2$$
We can then apply the Euler-Lagrange equation to both θ₁ and θ₂:
$$\frac{\partial L}{\partial \theta_1} = \frac{d}{dt} \left(\frac{\partial L}{\partial \dot{\theta_1}}\right)
\implies \ddot\theta_1 + g_1(\theta_1, \theta_2)\ddot\theta_2 = f_1(\theta_1, \theta_2, \dot\theta_1, \dot\theta_2)$$
$$\frac{\partial L}{\partial \theta_2} = \frac{d}{dt} \left(\frac{\partial L}{\partial \dot{\theta_2}}\right)
\implies \ddot\theta_2 + g_2(\theta_1, \theta_2)\ddot\theta_1 = f_2(\theta_1, \theta_2, \dot\theta_1, \dot\theta_2)$$
where for readability I have defined functions:
$$ g_1 = \frac{m_2}{M}\frac{L_2}{L_1}\cos\Delta\theta, \hspace{3mm}
f_1 = -\frac{m_2}{M}\frac{L_2}{L_1}{\dot\theta_2}^{2}\sin\Delta\theta - \frac{g}{L_1}\sin\theta_1 $$
$$ g_2 = \frac{L_1}{L_2}\cos\Delta\theta, \hspace{3mm}
f_2 = \frac{L_1}{L_2}{\dot\theta_1}^{2}\sin\Delta\theta - \frac{g}{L_2}\sin\theta_2$$
We can express this system of equations in matrix form:
$$
\boldsymbol{G}\ddot{\vec{\theta}} =
\begin{pmatrix}
1 & g_1\\
g_2 & 1
\end{pmatrix}
\begin{pmatrix}
\ddot\theta_1 \\
\ddot\theta_2
\end{pmatrix} =
\begin{pmatrix}
f_1 \\
f_2
\end{pmatrix}
$$
To solve we note that:
$$\det(\boldsymbol{G}) = 1 - g_1g_2 = 1 - \frac{m_2}{M}\cos^{2}\Delta\theta \ge 1 - \frac{m_2}{M} > 0$$
Thus, G is a non-singular matrix and has an inverse. We can apply G inverse to both sides to solve:
$$\ddot{\vec{\theta}} = \boldsymbol{G}^{-1}\vec{f} = \frac{1}{\det(\boldsymbol{G})}
\begin{pmatrix}
1 & -g_1 \\
-g_2 & 1
\end{pmatrix}
\begin{pmatrix}
f_1 \\
f_2
\end{pmatrix} =
\frac{1}{1-g_1g_2}
\begin{pmatrix}
f_1 - g_1f_2 \\
-g_2f_1 + f_2
\end{pmatrix} =
\begin{pmatrix}
h_1(\theta_1, \theta_2, \dot\theta_1, \dot\theta_2) \\
h_2(\theta_1, \theta_2, \dot\theta_1, \dot\theta_2)
\end{pmatrix}
$$
where again for readability I have defined functions:
$$ h_1 = \frac{f_1 - g_1f_2}{1-g_1g_2}$$
$$ h_2 = \frac{-g_2f_1 + f_2}{1-g_1g_2}$$
Finally, to get it into a form useful for numerical integration we turn it into
a set of first-order differential equations by defining ωᵢ = θᵢ':
$$\frac{d}{dt}
\begin{pmatrix}
\theta_1 \\
\theta_2 \\
\omega_1 \\
\omega_2
\end{pmatrix} =
\begin{pmatrix}
\omega_1 \\
\omega_2 \\
h_1 \\
h_2
\end{pmatrix}$$
This is a set of coupled non-linear ordinary differential equations and
cannot be solved analytically. However, we can use this form for
numerical integration. In my simulation I use the fourth-order Runge-Kutta method (RK4).
Numerical Integration Method
This simulation uses the fourth-order Runge-Kutta (RK4) method for numerically solving
the system of first-order non-linear ordinary differential equations presented above.
Since the double pendulum is an autonomous system the initial value problem can be specified as follows:
$$ \frac{d\vec{y}}{dt} = \vec{f}(\vec{y}), \hspace{3mm} \vec{y}(t_0) = \vec{y_0}$$
where:
$$ \vec{y} = \begin{pmatrix}
\theta_1 \\
\theta_2 \\
\omega_1 \\
\omega_2
\end{pmatrix}, \hspace{3mm} \vec{f} =
\begin{pmatrix}
\omega_1 \\
\omega_2 \\
h_1 \\
h_2
\end{pmatrix}
$$
The Runge-Kutta method takes this IVP and a time step, h, and
determines the next state based on the previous state and a successive weighted average of slopes:
$$\vec{y}_{n+1} = \vec{y}_n + \frac{h}{6}(\vec{k_1} + 2\vec{k_2} + 2\vec{k_3} + \vec{k_4})$$
$$t_{n+1} = t_n + h$$
where the slopes are:
$$\vec{k_1} = \vec{f}(\vec{y_n})$$
$$\vec{k_2} = \vec{f}(\vec{y_n} + \frac{h}{2}\vec{k_1})$$
$$\vec{k_3} = \vec{f}(\vec{y_n} + \frac{h}{2}\vec{k_2})$$
$$\vec{k_4} = \vec{f}(\vec{y_n} + h\vec{k_3})$$