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:

peng_drawing

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! :D

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!

q1_cartdisp

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.

q1_energy

The energy graph shows that kinetic energy (KE) added to potential energy (PE) results in a flat total energy (TE). Energy is conserved, yay!

q1_m1_dva

The position/velocity/acceleration plot for M1 is quite turbulent, because M1 has more forces acting on it than M2.

q1_m2_dva

M2’s motion is as expected.

q1_tension

Finally, we see that the tension forces on M1 and M2 add to 0. Because Newton’s third law!

Whew, that was quite a ride. Cool, innit?

-Sophie

Leave a Reply

Your email address will not be published. Required fields are marked *