Symbolic Math Toolbox

Double Pendulum Modeling

Introduction

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

image

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

math
math

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

math
math
math
math

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

math
math
math
math

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

math
math
math
math

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.

image

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;

math
math

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.

image

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;

math
math

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)

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

math
Click on image to see enlarged view.
math
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):

math
Initial Conditions:
math

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

MuPAD graphics

製品評価版の入手または製品の購入