The double compound pendulum

Posted by: christian on 15 Apr 2020

(6 comments)

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.

double compound pendulum animation

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)
    ax.set_aspect('equal', adjustable='box')
    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)
ax = fig.add_subplot(111)

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.

Geometry for the derivation of the dynamical equations for the double compound pendulum

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

Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

Peter Thorsteinson 4 years, 7 months ago

I get "Maximum energy drift exceeded"

Link | Reply
Current rating: 5

christian 4 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.

Link | Reply
Current rating: 5

Jonathan 1 year, 6 months ago

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

Link | Reply
Current rating: 5

Christian Hill 1 year, 6 months ago

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

Link | Reply
Current rating: 5

Donato Paul 9 months, 2 weeks ago

Really nice article about a double complex pendulum, it was really helpful (helped a lot in verifying my own Lagrangian)!

I just have one doubt about the moments of inertia used to get the rods' rotational energy: Aren't the axes of rotation at the endpoint of the rods? So, shouldn't the moment be equal to M^2*L/3, and not M^2*L/12 (when the axis of rotation is at the center of the rod)?

Link | Reply
Current rating: 5

christian 9 months, 2 weeks ago

Thank you! I think the issue you are referring to is discussed here: https://physics.stackexchange.com/questions/62618/double-compound-pendulum-why-use-inertia-about-the-center-of-mass-for-bottom-pe and in the link to the chat section of the Wikipedia article – hope that helps.
Cheers, Christian

Link | Reply
Current rating: 5

New Comment

required

required (not published)

optional

required