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)

We begin by defining symbols for θ_{1} and θ_{2}.

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

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

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

We define displacement expressions for m_{1} and m_{2 } 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));

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);

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);

a_y_1 := t --> v_y_1'(t);

a_x_2 := t --> v_x_2'(t);

a_y_2 := t --> v_y_2'(t);

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 (a_{x1} and a_{x2}) 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;

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 m_{2}. Our free body diagram shows that there are 2 forces acting on m_{2} ; 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;

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

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

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

eq5 := t --> subs(eq3,op(Tension[1],1..2));

eq6 := t --> subs(eq4,op(Tension[1],1..2));

eq6 := t --> subs(eq4,op(Tension[1],1..2));

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:

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

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

);

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

);