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:
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θ2Note: 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=˙y1−l2˙θ2sinθ2Our Lagrangian is of the form T−V 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=−m1gy1−m2gy2=−(m1+m2)gl1cosθ1−m2gl2cosθ2Now we can compute our Lagrangian L:
L=T−V=12l21˙θ21(m1+m2)+12m2(2l1˙θ1l2˙θ2cos(θ1−θ2)+l22˙θ22)+(m1+m2)gl1cosθ1+m2gl2cosθ2Time 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θ1And 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¨θ2∂L∂θ2=m2l1l2˙θ1˙θ2sin(θ1−θ2)−m2gl2sinθ20=∂L∂θ2−˙Pθ2=m2l2¨θ2+m2l1(¨θ1cos(θ1θ2)−˙θ21sin(θ1−θ2))+m2gsinθ2Solving for ¨θ1 and ¨θ2 we get:
¨θ1=gm1sin(θ1)+gm22sin(θ1)+gm22sin(θ1−2θ2)+l1m22˙θ21sin(2θ1−2θ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θ1−2θ2)2l2(m1−m2cos2(θ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(θ1−2θ2)+l1m22ω21sin(2θ1−2θ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θ1−2θ2)2l2(m1−m2cos2(θ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)=0Here’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.