# The double compound pendulum

Following on from this post about the simple double pendulum, (two bobs connected by light, rigid rods), this post animates the double compound pendulum (also called a double complex or physical pendulum): two rods connected to each other, with their mass distributed along their length. The analysis on Wikipedia provides the dynamical equations for the case of equal-mass and equal-length rods. Here, the more general case of rods with lengths $l_1$ and $l_2$ and masses $m_1$ and $m_2$ is considered. The code below requires SciPy 1.4+.

import sys
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Pendulum rod lengths (m) and masses (kg).
L1, L2 = 1, 1.5
m1, m2 = 1, 2
# The gravitational acceleration (m.s-2).
g = 9.81
# Check that the integration conserves total energy to within this absolute
# tolerance.
EDRIFT = 0.05

def deriv(t, y, L1, L2, m1, m2):
"""Return the derivatives of y = theta1, p_theta1, theta2, p_theta2.

These are the generalized coordinates (here, angles) and generalized
momenta for the two rigid rods.

"""

theta1, p_theta1, theta2, p_theta2 = y

M = m1 + 3*m2
c, s = np.cos(theta1 - theta2), np.sin(theta1 - theta2)
Lr = L1 / L2
den = 4 * M - 9 * m2 * c**2

theta1dot = 6 / L1**2 * (2*p_theta1 - 3 * Lr * c * p_theta2) / den
theta2dot = 6 / m2 / L2**2 * (
(2 * p_theta2 * M - 3 * m2 / Lr * c * p_theta1) / den)
term = m2 * L1 * L2 / 2 * theta1dot * theta2dot * s
p_theta1dot = -term - (m1/2 + m2) * g * L1 * np.sin(theta1)
p_theta2dot = term - m2/2 * g * L2 * np.sin(theta2)

return theta1dot, p_theta1dot, theta2dot, p_theta2dot

def calc_H(y, L1, L2, m1, m2):
"""Calculate the Hamiltonian at y = theta1, p_theta1, theta2, p_theta2."""
theta1, p_theta1, theta2, p_theta2 = y

theta1dot, p_theta1dot, theta2dot, p_theta2dot = deriv(None, y, L1, L2,
m1, m2)
# The Lagrangian
c = np.cos(theta1 - theta2)
L = ( m1 * (L1 * theta1dot)**2 / 6 + m2 * (L2 * theta2dot)**2 / 6
+ m2 / 2 * ((L1 * theta1dot)**2 + L1*L2*theta1dot*theta2dot * c)
+ g * L1 * np.cos(theta1) * (m1 / 2 + m2)
+ g * L2 * np.cos(theta2) * m2 / 2
)

# The Hamiltonian
H = p_theta1 * theta1dot + p_theta2 * theta2dot - L
return H

# Maximum time, time point spacings and the time grid (all in s).
tmax, dt = 10, 0.01
t = np.arange(0, tmax+dt, dt)

# Initial conditions:
# angles: theta1, theta2 and generalized momenta: p_theta1, p_theta2
theta1_0, theta2_0 = np.pi, np.pi/2
p_theta1_0, p_theta2_0 = 0, 0
y0 = np.array([theta1_0, p_theta1_0, theta2_0, p_theta2_0])
# Could call calc_H, but since the initial p_thetai are zero, and the
# Hamiltonian is conserved (since the Langrangian has no explicit time-
# dependence, H0 is just the potential energy of the initial configuration:
H0 = -g * (L1 * np.cos(theta1_0) * (m1 / 2 + m2) +
L2 * np.cos(theta2_0) * m2 / 2)

# Do the numerical integration of the equations of motion.
y = solve_ivp(deriv, (0, tmax), y0, method='Radau', dense_output=True,
args=(L1, L2, m1, m2))

# Check that the Hamiltonian didn't drift too much.
H = calc_H(y.y, L1, L2, m1, m2)
if any(abs(H-H0) > EDRIFT):
print('Maximum energy drift exceeded')

# Unpack dynamical variables as a function of time.
theta1, p_theta1, theta2, p_theta2 = y.sol(t)

# Convert to Cartesian coordinates of the two rods.
x1 = L1 * np.sin(theta1)
y1 = -L1 * np.cos(theta1)
x2 = x1 + L2 * np.sin(theta2)
y2 = y1 - L2 * np.cos(theta2)

# Plot a trail of the m2 bob's position for the last trail_secs seconds.
trail_secs = 1
# This corresponds to max_trail time points.
max_trail = int(trail_secs / dt)

def make_plot(i):
"""
Plot and save an image of the double pendulum configuration for time
point i.
"""

# The pendulum rods (thick, black).
ax.plot([0, x1[i], x2[i]], [0, y1[i], y2[i]], lw=8, c='k')

# The trail will be divided into ns segments and plotted as a fading line.
ns = 20
s = max_trail // ns

for j in range(ns):
imin = i - (ns-j)*s
if imin < 0:
continue
imax = imin + s + 1
# The fading looks better if we square the fractional length along the
# trail.
alpha = (j/ns)**2
# Two trails, initiating at the centres of mass of the two rods.
ax.plot(x1[imin:imax]/2, y1[imin:imax]/2, c='b', solid_capstyle='butt',
lw=2, alpha=alpha)
ax.plot((x1[imin:imax]+x2[imin:imax])/2,
(y1[imin:imax]+y2[imin:imax])/2,
c='r', solid_capstyle='butt', lw=2, alpha=alpha)

# Centre the image on the fixed anchor point, and ensure the axes are equal
ax.set_xlim(-L1-L2, L1+L2)
ax.set_ylim(-L1-L2, L1+L2)
plt.axis('off')
plt.savefig('frames/_img{:04d}.png'.format(i//di), dpi=72)
plt.cla()

# Make an image every di time points, corresponding to a frame rate of fps
# frames per second.
# Frame rate, s-1
fps = 10
di = int(1/fps/dt)
fig = plt.figure(figsize=(8.3333, 6.25), dpi=72)

for i in range(0, t.size, di):
print(i // di, '/', t.size // di)
make_plot(i)


### Derivation of the dynamical equations

This derivation follows the Hamiltonian formulation of the dynamics of the double compound pendulum and is slightly more general than that given in this Wikipedia page.

The pendulum rods are of lengths $l_1$ and $l_2$ and have masses $m_1$ and $m_2$ uniformly distributed along their lengths. The coordinate system used is illustrated below. In terms of the angles $\theta_1$ and $\theta_2$, the centres of mass of the rods are at the coordinates:

\begin{align*} (x_1, y_1) &= \textstyle(\frac{1}{2}l_1\sin\theta_1, \; -\frac{1}{2}l_1\cos\theta_1),\\ (x_2, y_2) &= \textstyle(l_1\sin\theta_1 + \frac{1}{2}l_2\sin\theta_2, \; -l_1\cos\theta_1 - \frac{1}{2}l_2\cos\theta_2). \end{align*}

Their first derivatives with respect to time are therefore:

\begin{align*} (\dot{x}_1, \dot{y}_1) &= \textstyle(\frac{1}{2}l_1\dot{\theta}_1\cos\theta_1, \; \frac{1}{2}l_1\dot{\theta}_1\sin\theta_1),\\ (\dot{x}_2, \dot{y}_2) &= \textstyle(l_1\dot{\theta}_1\cos\theta_1 + \frac{1}{2}l_2\dot{\theta}_2\cos\theta_2, \; l_1\dot{\theta}_1\sin\theta_1 + \frac{1}{2}l_2\dot{\theta}_2\sin\theta_2). \end{align*}

The total kinetic energy of each rod is most conveniently thought of as being composed of a translational part, $\frac{1}{2}m_i(\dot{x}_i^2 + \dot{y}_i^2)$ and a rotational part, $\frac{1}{2}I_i\dot{\theta}_i^2$, where the moment of inertia of a thin, rigid rod of length $l_i$ and mass $m_i$ is $I_i = \frac{1}{12}m_il_i^2$. Therefore,

\begin{align*} T &= \textstyle \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) + \frac{1}{2}I_1\dot{\theta}_1^2 + \frac{1}{2}I_2\dot{\theta}_2^2\\ &= \textstyle \frac{1}{6}m_1l_1^2\dot{\theta}_1^2 + \frac{1}{6}m_2l_2^2\dot{\theta}_2^2 + \frac{1}{2}m_2[l_1^2\dot{\theta}_1^2 + l_1l_2\dot{\theta}_1\dot{\theta}_2\cos(\theta_1 - \theta_2)]. \end{align*}

The potential energy is $$V = m_1gy_1 + m_2gy_2 =\textstyle -gl_1(\frac{1}{2}m_1+m_2)\cos\theta_1 - \frac{1}{2}m_2gl_2\cos\theta_2,$$ and from these two quantities, the Lagrangian, $\mathcal{L} = T - V$ can be calculated.

The generalized momenta, $p_i$ are defined through the relation $$p_i = \frac{\partial\mathcal{L}}{\partial\dot{\theta}_i},$$ and so:

\begin{align*} p_1 &= \textstyle \frac{1}{3}m_1l_1^2\dot{\theta}_1 + m_2l_1^2\dot{\theta}_1 + \frac{1}{2}m_2l_1l_2\dot{\theta}_2\cos(\theta_1 - \theta_2),\\ p_2 &= \textstyle \frac{1}{3}m_2l_2^2\dot{\theta}_2 + \frac{1}{2}m_2l_1l_2\dot{\theta}_1\cos(\theta_1 - \theta_2). \end{align*}

With a bit of effort, these equations can be inverted to give the following expressions for the angular velocities:

\begin{align*} \dot{\theta}_1 &= \frac{6}{l_1^2}\left[ \frac{2p_1 - 3 \frac{l_1}{l_2}\cos(\theta_1 - \theta_2)p_2}{4(m_1+3m_2)-9m_2\cos^2(\theta_1 - \theta_2)}\right],\\ \dot{\theta}_2 &= \frac{6}{m_2l_2^2}\left[ \frac{2p_2(m_1+3m_2) - 3m_2 \frac{l_2}{l_1}\cos(\theta_1 - \theta_2)p_1}{4(m_1+3m_2)-9m_2\cos^2(\theta_1 - \theta_2)}\right]. \end{align*}

Now, the Euler-Lagrange equations,

\begin{align*} \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial\mathcal{L}}{\partial \dot{q}_i}\right) - \frac{\partial \mathcal{L}}{\partial q_i} = 0 \quad \mathrm{for}\;q_i = \theta_1, \theta_2, \end{align*}

imply that

$$\dot{p}_i = \frac{\partial\mathcal{L}}{\partial\theta_i},$$

and so:

\begin{align*} \dot{p}_1 &= \textstyle -\frac{1}{2}m_2l_1l_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1 - \theta_2) - (\frac{1}{2}m_1+m_2)gl_1\sin\theta_1,\\ \dot{p}_2 &= \textstyle \frac{1}{2}m_2l_1l_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1 - \theta_2) - \frac{1}{2}m_2gl_2\sin\theta_2.\\ \end{align*}

The equations giving $\dot{\theta}_1$, $\dot{\theta}_2$, $\dot{p}_1$ and $\dot{p}_2$ are all we need to model the time-evolution of the pendulum system. Since the Lagrangian does not have an explicit time-dependence, the Hamiltonian $H = T + V$ is conserved as the total energy, so the numerical integration of these equations can be checked to ensure that this quantity is constant.

Current rating: 5.0 #### Peter Thorsteinson 3 years, 7 months ago

I get "Maximum energy drift exceeded"

Current rating: 5 #### christian 3 years, 7 months ago

What were your initial conditions? It's possibly the integrator getting lost or a bug in the algorithm... but I did check the maths and thought it was OK.

Currently unrated #### Jonathan 7 months ago

In the position vector of your first pendulum rod, the y-component should be -l1/2*cos(phi1) instead of l2 right?

Currently unrated #### Christian Hill 6 months, 3 weeks ago

Yup – thanks. I corrected this typo.
Cheers, Christian