商店 
Symbolic Math Toolbox 5.5
Modeling the Motion of a Double Pendulum
Introduction
We will model the dynamic behavior of a double pendulum.
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 Masses
We begin by defining symbols for θ1 and θ2.
th2 := Symbol::subScript(Symbol("theta"),2);
We define displacement expressions for m1 and m2 and calculate velocity by differentiating displacement with respect to time.
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_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_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 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.
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.
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.
Click on image to see enlarged view.
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_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.
// 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.
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
);
