# Symbolic Math Toolbox

## Double Pendulum Modeling

### Introduction

This example shows how to model the motion of a double pendulum using symbolic computation.

where:

x - horizontal position of pendulum mass

y - vertical position of pendulum mass

θ - angle of pendulum (0 = vertical downwards, counter-clockwise is positive)

L - length of rod (constant)

### Define Displacement, Velocity, and Acceleration of Double Pendulum Masses

We begin by defining symbols for θ1 and θ2.

th1 := Symbol::subScript(Symbol("theta"),1);
th2 := Symbol::subScript(Symbol("theta"),2);

We define displacement expressions for m1 and m2 and calculate velocity by differentiating displacement with respect to time.

x_1 := t --> L_1 * sin(th1(t));
y_1 := t --> - L_1 * cos(th1(t));
x_2 := t --> x_1(t) + L_2 * sin(th2(t));
y_2 := t --> y_1(t) - L_2 * cos(th2(t));

v_x_1 := t --> x_1'(t);
v_y_1 := t --> y_1'(t);
v_x_2 := t --> x_2'(t);
v_y_2 := t --> y_2'(t);

We calculate acceleration by differentiating velocity with respect to time.

a_x_1 := t --> v_x_1'(t);
a_y_1 := t --> v_y_1'(t);
a_x_2 := t --> v_x_2'(t);
a_y_2 := t --> v_y_2'(t);

### Evaluate Forces on Double Pendulum Masses

We will now evaluate the forces acting on the masses. We start by constructing a free body diagram of the forces on m1.

Balancing the horizontal and vertical force components results in 2 equations.  Note that we use the acceleration expressions that we calculated previously (ax1 and ax2) in our equations.

eq1 := m_1*a_x_1(t) = − T_1*sin(th1(t)) + T_2*sin(th2(t));
eq2 := m_1*a_y_1(t) = T_1*cos(th1(t)) − T_2*cos(th2(t)) − m_1*g;

We now evaluate the forces acting on m2.  Our free body diagram shows that there are 2 forces acting on m2 ; tension from rod 2, and the force from the weight of the mass itself.

Balancing the horizontal and vertical force components results in 2 equations.

eq3 := m_2*a_x_2(t) = − T_2*sin(th2(t));
eq4 := m_2*a_y_2(t) = T_2*cos(th2(t)) − m_2*g;

### Reduce System Equations

Our force evaluation produced a set of 4 equations with 4 unknowns ([θ1 θ2 T1 T2]). We will reduce our system to 2 equations with 2 unknowns by solving eq1 and eq2 for T1 and T2, and substituting the results into eq3 and eq4.

Tension := solve([eq1,eq2],[T_1,T_2],IgnoreSpecialCases)

Click on image to see enlarged view.
eq5 := t --> subs(eq3,op(Tension[1],1..2));
eq6 := t --> subs(eq4,op(Tension[1],1..2));

Click on image to see enlarged view.

Click on image to see enlarged view.

### Solve System Equations and Animate Pendulum Motion

Before solving the system equations, we define the masses and rod lengths.

L_1 := 1:
L_2 := 1.5:
m_1 := 1:
m_2 := 2:
g := 98/10:

We specify initial conditions for mass position and angular velocity, and solve the final system equations (eq5, eq6) numerically using the numeric::odesolve2 function.

fields := [th1(t),th1'(t),th2(t),th2'(t)];

// Initial conditions
print(Unquoted, "Initial Conditions:");
Y0 := [PI/6, 0, PI/4, 0];

IVP := {eq5(t),
eq6(t),
th1(0)  = Y0[1],
th1'(0) = Y0[2],
th2(0)  = Y0[3],
th2'(0) = Y0[4]
}:
odesolve2args := numeric::ode2vectorfield(IVP, fields):
Ynum := numeric::odesolve2(odesolve2args):

Initial Conditions:

We animate the motion of the double pendulum.  Note that this animation will not run in a Web browser, however you can download this example from MATLAB Central (http://www.mathworks.com/matlabcentral/) and run it yourself.

dt := 0.05:
imax := 60:
plot(
plot::Line2d([0, 0], [L_1 * sin(Ynum(t)[1]), -L_1 * cos(Ynum(t)[1])], VisibleFromTo = t..t + 0.99*dt,
LineColor = RGB::Red) $t in [i*dt$ i = 0..imax],
plot::Point2d([L_1 * sin(Ynum(t)[1]), -L_1 * cos(Ynum(t)[1])], VisibleFromTo = t..t + 0.99*dt, PointSize = m_1*4*unit::mm,
Color = RGB::Red) $t in [i*dt$ i = 0..imax],
plot::Line2d([L_1 * sin(Ynum(t)[1]), -L_1 * cos(Ynum(t)[1])], [L_1 * sin(Ynum(t)[1]) + L_2 * sin(Ynum(t)[3]), -L_1 * cos(Ynum(t)[1]) - L_2 * cos(Ynum(t)[3])], VisibleFromTo = t..t + 0.99*dt,
LineColor = RGB::Green) $t in [i*dt$ i = 0..imax],
plot::Point2d([L_1 * sin(Ynum(t)[1]) + L_2 * sin(Ynum(t)[3]), -L_1 * cos(Ynum(t)[1]) - L_2 * cos(Ynum(t)[3])], VisibleFromTo = t..t + 0.99*dt, PointSize = m_2*4*unit::mm,
Color = RGB::Green) $t in [i*dt$ i = 0..imax],
Header = "Double Pendulum",
AnimationStyle = Loop
);