Equations of motion for a planar simple double pendulum

To provide some background information for my N-link pendulum project, I’ve broken the methodology for solving the equations of motion (EOM) for a simple double pendulum into a separate post. I’ll explain to the best of my abilities!

The pendulum that we are solving for is composed of two masses, m1 and m2, suspended from each other by strings of length l1 and l2. Their positions relative to 0 is represented by theta1 and theta2. The entire pendulum is supported by point 0, which is a frictionless, massless point. The whole thing is a pain to create digitally, so I drew the setup below:

This will be a pretty long post, so scroll to the bottom for the cool graphs!

We’ll find the equations of motion in Polar coordinates, since it means that we only need two equations instead of four. (The two being $\ddot{\theta_1}, \ddot{\theta_2}$ and the four being ax1, ay1, ax2, ay2.)

By doing basic trig, we can find the EOM of the masses using time derivatives of the unit vectors. In this case, $\dot{\theta_1} = \dfrac{d\theta_1}{dt}$ and $\ddot{\theta_1} = \dfrac{d^2\theta_1}{dt^2}$ .

Position:
$x_1 = L_1sin(\theta_1)$
$y_1 = -L_1cos(\theta_1)$

$x_2 = x_1 + L_2sin(\theta_2)$
$y_2 = y_1 - L_2cos(\theta_2)$

Velocity:
$vx_1 = L_1\dot{\theta_1}cos(\theta_1)$
$vy_1 = L_1\dot{\theta_1}sin(\theta_1)$

$vx_2 = vx_1 + L_2\dot{\theta_2}cos(\theta_2)$
$vy_2 = vy_1 + L_2\dot{\theta_2}sin(\theta_2)$

Acceleration:
$ax_1 = L_1\ddot{\theta_1}cos(\theta_1) - L_1\dot{\theta_1}^2sin(\theta_1)$
$ay_1 = L_1\ddot{\theta_1}sin(\theta_1) + L_1\dot{\theta_1}^2cos(\theta_1)$

$ax_2 = ax_1 + L_2\ddot{\theta_2}cos(\theta_2) - L_2\dot{\theta_2}^2sin(\theta_2)$
$ay_2 = ay_1 + L_2\ddot{\theta_2}sin(\theta_2) + L_2\dot{\theta_2}^2cos(\theta_2)$

We know that $F = ma$ , so to properly simulate this system we must also account for the tension and gravity forces.

The sum of forces is:
$-T_1sin(\theta_1) + T_2sin(\theta_2) = m_1 * ax_1$
$T_1cos(\theta_1) - T_2cos(\theta_2) - m_1g = m_1 * ay_1$

$-T_2sin(\theta_2) = m_2 * ax_2$
$T_2cos(\theta_2) - m_2g = m_2 * ay_2$

Note at ax1, ay1, etc. are the same as the acceleration equations we previously found.
Because tension is a reactionary force, we’ll solve for T1 and T2 in order to isolate equations for $\ddot{\theta_1}, \ddot{\theta_2}$ .

Solving for T1:
$-T_1sin(\theta_1) + T_2sin(\theta_2) = m_1 * ax_1$
$-T_1sin(\theta_1) = m_1 * ax_1 + m_2 * ax_2$
$T_1 = -\dfrac{m_1 * ax_1 + m_2 * ax_2}{sin(\theta_1)}$

$T_1cos(\theta_1) - T_2cos(\theta_2) - m_1g = m_1 * ay_1$
$T_1cos(\theta_1) = m_1 * ay_1 + m_2 * ay_2 + m_1g + m_2g$
$T_1 = \dfrac{m_1 * ay_1 + m_2 * ay_2 + m_1g + m_2g}{cos(\theta_1)}$

Setting the two equations equal to each other:
$-cos(\theta_1)[m_1 * ax_1 + m_2 * ax_2] = sin(\theta_1)[m_1 * ay_1 + m_2 * ay_2 + m_1g + m_2g]$

Now for T2:
$-T_2sin(\theta_2) = m_2 * ax_2$
$T_2 = -\dfrac{m_2 * ax_2}{\theta_2}$

$T_2cos(\theta_2) - m_2g = m_2 * ay_2$
$T_2 = \dfrac{m_2 * ay_2 + m_2g}{cos(\theta_2)}$

Setting the two equal:
$cos(\theta_2)m_2 * ax_2 = sin(\theta_2)[m_2 * ay_2 + m_2g]$

Now we have two equations and two unknowns! We can solve for $\ddot{\theta_1}, \ddot{\theta_2}$ in terms of $\theta_1, \dot{\theta_1}, \theta_2, \dot{\theta_2}$ . This was a bit of a behemoth venture and it took me 14 pages….

I did all the wrong math so you don’t have to! 😀

UNTIL I REALIZED that I had dropped a negative sign. Throw it into Matlab using the syms function and have it spit out the correct equations. Which, if you were curious, look like this:

t1dd = -(L1*M2*sin(2*t1 - 2*t2)*t1d^2 + 2*L2*M2*sin(t1 - t2)*t2d^2 + M2*g*sin(t1 - 2*t2) + 2*M1*g*sin(t1) + M2*g*sin(t1))/(L1*(2*M1 + M2 - M2*cos(2*t1 - 2*t2)));

t2dd = (M1*g*sin(2*t1 - t2) + M2*g*sin(2*t1 - t2) - M1*g*sin(t2) - M2*g*sin(t2) + L2*M2*t2d^2*sin(2*t1 - 2*t2) + 2*L1*M1*t1d^2*sin(t1 - t2) + 2*L1*M2*t1d^2*sin(t1 - t2))/(L2*(2*M1 + M2 - M2*cos(2*t1 - 2*t2)));

(terrifying).

Using all this information, you can put the equations into Matlab’s ODE45 to plot the motion of the simple pendulum!

The cartesian displacement of the two masses. Note how M1 is circular (as expected), while M2 moves in a bowl like path because its trajectory is dependent on M1.