double pendulum motion

Alex Vig bio photo By Alex Vig

intro

This is an implementation of the double pendulum example given in Leonard Susskind’s lecture on classical mechanics. The goal of this post is to explain the derivation and implementation of the equations of motion for coupled pendulums.

The accompanying Python implementation is available on GitHub.

I’ll start with the drawing given in Professor Susskind’s lecture with one exception. The angle α will be called θ2 and it will be the angle with respect to the rod securing the apparatus (same relationship as θ1) rather than the angle with respect to θ1.

So here’s our system:

System Image courtesy of JabberWok.

derivation

Using the image above and some basic trigonometry we can derive the following equations for the x, y positions of the pendulums:

x1=l1sinθ1y1=l1cosθ1x2=x1+l2sinθ2y2=y1+l2cosθ2

Note: In this inertial reference frame the positive y direction is down.

Differentiating with respect to time we get:

˙x1=l1˙θ1cost1˙y1=l1˙θ1sint1˙x2=˙x1+l2˙θ2cosθ2˙y2=˙y1l2˙θ2sinθ2

Our Lagrangian is of the form TV where T is the kinetic energy and V is the potential energy.

Let’s calculate the kinetic energy T:

T=12m1(˙x21+˙y21)+12m2(˙x22+˙y22)=12m1l21˙θ21+12m2(l21˙θ21+2l1˙θ1l2˙θ2(cosθ1cosθ2+sinθ1sinθ2)+l22˙θ22)=12l21˙θ21(m1+m2)+12m2(2l1˙θ1l2˙θ2cos(θ1θ2)+l22˙θ22)

And now the potential energy V:

V=m1gy1m2gy2=(m1+m2)gl1cosθ1m2gl2cosθ2

Now we can compute our Lagrangian L:

L=TV=12l21˙θ21(m1+m2)+12m2(2l1˙θ1l2˙θ2cos(θ1θ2)+l22˙θ22)+(m1+m2)gl1cosθ1+m2gl2cosθ2

Time to get the terms for the Euler-Lagrange equation for the coordinate θ1:

Pθ1=l21˙θ1(m1+m2)+m2l1l2˙θ2cos(θ1θ2)˙Pθ1=l21¨θ1(m1+m2)+m2l1l2(¨θ2cos(θ1θ2)˙θ2sin(θ1θ2)(˙θ1˙θ2))Lθ1=m2l1td1l2td2sin(θ1θ2)(m1+m2)gl1sinθ10=Lθ1˙Pθ1=l1tdd1(m1+m2)+m2l2(tdd2cos(θ1θ2)+td22sin(θ1θ2))+(m1+m2)gsinθ1

And now the terms for the Euler-Lagrange equation for the coordinate θ2:

Pθ2=m2l1l2˙θ1cos(θ1θ2)+m2l22˙θ2˙Pθ2=m2l1l2(¨θ1cos(θ1θ2)˙θ1sin(θ1θ2)(˙θ1˙θ2))+m2l22¨θ2Lθ2=m2l1l2˙θ1˙θ2sin(θ1θ2)m2gl2sinθ20=Lθ2˙Pθ2=m2l2¨θ2+m2l1(¨θ1cos(θ1θ2)˙θ21sin(θ1θ2))+m2gsinθ2

Solving for ¨θ1 and ¨θ2 we get:

¨θ1=gm1sin(θ1)+gm22sin(θ1)+gm22sin(θ12θ2)+l1m22˙θ21sin(2θ12θ2)+l2m2˙θ22sin(θ1θ2)l1(m1+m2cos2(θ1θ2)m2)¨θ2=gm1sin(θ2)+gm1sin(2θ1θ2)gm2sin(θ2)+gm2sin(2θ1θ2)+2l1m1˙θ21sin(θ1θ2)+2l1m2˙θ21sin(θ1θ2)+l2m2˙θ22sin(2θ12θ2)2l2(m1m2cos2(θ1θ2)+m2)

Note: I didn’t do some crazy mental variable substitutions/arithmetic. I got these expressions using sympy:

import sympy

s = 'l_1 theta_dd_1 m_1 m_2 l_2 theta_dd_2 theta_1 theta_2 theta_d_2 theta_d_1 g'
l_1, theta_dd_1, m_1, m_2, l_2, theta_dd_2, theta_1, theta_2, theta_d_2, theta_d_1, g = sympy.symbols(s)
expr1 = l_1 * theta_dd_1 * (m_1 + m_2) + m_2 * l_2 * (theta_dd_2 * sympy.cos(theta_1 - theta_2) + theta_d_2 ** 2 * sympy.sin(theta_1 - theta_2)) + (m_1 + m_2) * g * sympy.sin(theta_1)
expr2 = m_2 * l_2 * theta_dd_2 + m_2 * l_1 * (theta_dd_1 * sympy.cos(theta_1 - theta_2) - theta_d_1 ** 2 * sympy.sin(theta_1 - theta_2)) + m_2 * g * sympy.sin(theta_2)
theta_dd_1_expr = sympy.solveset(expr1, theta_dd_1).args[0]
theta_dd_2_eval = sympy.solveset(expr2.subs(theta_dd_1, theta_dd_1_expr), theta_dd_2).args[0]
theta_dd_2_expr = sympy.solveset(expr2, theta_dd_2).args[0]
theta_dd_1_eval = sympy.solveset(expr1.subs(theta_dd_2, theta_dd_2_expr), theta_dd_1).args[0]

print(sympy.latex(sympy.simplify(theta_dd_1_eval)))
print(sympy.latex(sympy.simplify(theta_dd_2_eval)))

solution

Now that we have coupled differential equations we can solve for motion. We will need to solve numerically and by reframing this system of second order ODEs as four first order differential equations we will be able to do so:

ω1=˙θ1ω2=˙θ2˙ω1=gm1sin(θ1)+gm22sin(θ1)+gm22sin(θ12θ2)+l1m22ω21sin(2θ12θ2)+l2m2ω22sin(θ1θ2)l1(m1+m2cos2(θ1θ2)m2)˙ω2=gm1sin(θ2)+gm1sin(2θ1θ2)gm2sin(θ2)+gm2sin(2θ1θ2)+2l1m1ω21sin(θ1θ2)+2l1m2ω21sin(θ1θ2)+l2m2ω22sin(2θ12θ2)2l2(m1m2cos2(θ1θ2)+m2)

We’ll use the fourth order Runge-Kutta method as implemented in the Scipy ODE solver:

def double_pend(
    y: list,
    t: float,
    m_1: float,
    l_1: float,
    m_2: float,
    l_2: float,
    g: float
) -> list:
    """
    Calculate one step for coupled pendulums.

    Args:
        y: (theta_1, theta_d_1, theta_2, theta_d_2) at this timestep.
        t: Time (unused).
        m_1: Mass 1.
        l_1: Length of pendulum 1.
        m_2: Mass 2.
        l_2: Length of pendulum 2.
        g: Gravity constant.

    Returns:
        Derivative estimates for each element of y.
    """

    theta_1, theta_d_1, theta_2, theta_d_2 = y

    theta_dd_1 = (g*m_1*sin(theta_1) + g*m_2*sin(theta_1)/2 + g*m_2*sin(theta_1 - 2*theta_2)/2 + l_1*m_2*theta_d_1**2*sin(2*theta_1 - 2*theta_2)/2 + l_2*m_2*theta_d_2**2*sin(theta_1 - theta_2))/(l_1*(-m_1 + m_2*cos(theta_1 - theta_2)**2 - m_2))
    theta_dd_2 = (-g*m_1*sin(theta_2) + g*m_1*sin(2*theta_1 - theta_2) - g*m_2*sin(theta_2) + g*m_2*sin(2*theta_1 - theta_2) + 2*l_1*m_1*theta_d_1**2*sin(theta_1 - theta_2) + 2*l_1*m_2*theta_d_1**2*sin(theta_1 - theta_2) + l_2*m_2*theta_d_2**2*sin(2*theta_1 - 2*theta_2))/(2*l_2*(m_1 - m_2*cos(theta_1 - theta_2)**2 + m_2))

    dydt = [theta_d_1, theta_dd_1, theta_d_2, theta_dd_2]

    return dydt


import functools
import numpy as np

system_params = {
    'm_1': 1,
    'l_1': 1,
    'm_2': 1,
    'l_2': 1,
    'g': 9.81
}

pend_func = functools.partial(double_pend, **system_params)
y0 = [np.pi / 2, 0, np.pi, 0]
t = np.linspace(0, 20, 1000)

theta_1, theta_d_1, theta_2, theta_d_2 = odeint(pend_func, y0, t).T

We can then convert these into x, y coordinates:

l_1, l_2 = system_params['l_1'], system_params['l_2']

x_1 = l_1 * np.sin(theta_1)
y_1 = l_1 * np.cos(theta_1)
x_2 = x_1 + l_2 * np.sin(theta_2)
y_2 = y_1 + l_2 * np.cos(theta_2)

Now we have an initial value problem so let’s pick some initial values:

θ1(0)=π2ω1(0)=0θ2(0)=πω2(0)=0

Here’s an animation of this system (plotted with matplotlib):

Feel free to checkout the code and play with the system by changing the parameters and initial values.

additional resources

  • All the code used to generate the double pendulum motions can be found on GitHub.
  • Leonard Susskind’s lecture that inspired this example and the course containing the lecture.
  • MyPhysicsLab solves the ODEs in javascript so you can change the initial conditions and see the corresponding motions in a browser.