Comparing the derivation of equation of motion for double pendulum by method of F=ma and by energy (Lagrangian) method

by Nasser M. Abbasi
May 11, 2011

The goal of this note is to show how to use a symbolic program to help solve a typical basic engineering problem that requires large amount of algebraic manipulation. By applying a CAS program such as Mathematica on such a problem, one can then concentrate more on the problem formulation itself by having the program do the long algebraic computation, which can be very time consuming if done by hand, and also can be error prone.
 
This also allows one to be able to experiment with different formulation, as one is now relieved from the time consuming part and is able to try different solutions and different idea, which can lead to better insight of the problem.

In this note, we solve the double pendulum problem both ways. One way is by using the energy method, and another way is by using F=ma approach. Both would require long time to do by hand due to the large amount of algebra involved. And at the end, both solutions are verified to be the same.

download   PDF  NOTEBOOK

Clear["Global`*"]
Needs["Notation`"]
Symbolize[(NotationTemplateTag[(Subscript[m, 1])])]
Symbolize[(NotationTemplateTag[(Subscript[m, 2])])]
Symbolize[(NotationTemplateTag[(Subscript[L, 1])])]
Symbolize[(NotationTemplateTag[(Subscript[L, 2])])]
rep[X_]:=X/.{r'[t]->(Overscript[r, •]),θ'[t]->(Overscript[θ, •]),θ[t]->θ,x'[t]->(Overscript[x, •]),r[r]->r,θ'[t]->(Overscript[θ, •]),r[t]->r,x[t]->x,ϕ'[t]->(Overscript[ϕ, •]),ϕ''[t]->(Overscript[ϕ, ••]),θ''[t]->(Overscript[θ, ••]),ϕ[t]->ϕ};

Simple pendulum

T=(1/2) m1*( L1* θ1'[t])^2;
V=m1 g L1 (1-Cos[ θ1[t]]);
L=T-V;
eq1=(D[D[L,θ1'[t]],t]-D[L,θ1[t]]==-c  θ1'[t])//FullSimplify

pars={L1->1,m1->1,g->9.81,c->0};
ic={θ1[0]==90 Degree  ,θ1'[0]==0};
θ1sol=θ1[t]/.First@NDSolve[Flatten[{eq1,ic}]/.pars,θ1[t],{t,0,20}];
θ1Speedsol=D[θ1sol,t]

Plot[θ1sol,{t,0,20}]

plot phase portrait

ParametricPlot[{θ1sol,θ1Speedsol},{t,0,20}]

another way to do the above

T=(1/2) m1*( L1* x1'[t])^2;
V=-m1 g L1 Cos[ x1[t]];
L=T-V;
eq1=(D[D[L,x1'[t]],t]-D[L,x1[t]]+c  x1'[t] )//FullSimplify
pars={L1->1,m1->1,g->9.81,c->1};
ic={x1[0]==90 Degree  ,x1'[0]==0};
sol=First@NDSolve[Flatten[{eq1==0,ic}]/.pars,x1[t],{t,0,20}]

der=D[sol,t]

Plot[x1[t]/.sol,{t,0,20},PlotRange->All]
ParametricPlot[Evaluate@{x1[t]/.sol,x1'[t]/.der},{t,0,20}]

sol

x1[t]/.sol/.t->3

D[x1[t]/.sol,t]/.t->4

θ1Speedsol/.t->3

der=D[sol,t]

Plot[x1'[t]/.der,{t,0,20},PlotRange->All]

system diagram

This diagram shows the velocity and acceleration components of the 2 masses in the double pendulum.

Lagrangian solution

Write down the kinetic energy T

T=(1/2)m1 (L1 θ1'[t])^2+
(1/2)m2 ( (L2 θ2'[t]Cos[θ2[t]]+L1 θ1'[t] Cos[θ1[t]])^2+
(L2 θ2'[t]Sin[θ2[t]]+L1 θ1'[t] Sin[θ1[t]])^2)

Simplify[%]

Row[{"T=",(1/2)(Subscript[m, 1]) ((Subscript[L, 1]) θ'[t])^2+FullSimplify[(1/2)(Subscript[m, 2]) ( ((Subscript[L, 2])(ϕ'[t])+(Subscript[L, 1]) θ'[t] Cos[ϕ[t]-θ[t]])^2+((Subscript[L, 1]) θ'[t] Sin[ϕ[t]-θ[t]])^2)]//rep}]

Write down the potential energy

V=m1 g( L1 (1-Cos[ θ1[t]]) )+m2 g(L1 (1-Cos[ θ1[t]])+ L2 (1-Cos[ θ2[t]]))

Vold=-L1 m1 g Cos[θ1[t]]-(L1 Cos[θ1[t]]+L2 Cos[θ2[t]])m2 g

Write down the Lagrangian L

L=T-V

Lold=T-Vold

L-Lold//Simplify

FullSimplify[L]

Find equation of motion for θ

eq1=D[D[L,θ1'[t]],t]-D[L,θ1[t]]+c θ1'[t]//Simplify

eq1old=D[D[Lold,θ1'[t]],t]-D[Lold,θ1[t]]+c θ1'[t]//Simplify

eq1-eq1old//FullSimplify

eq1=eq1/L1//FullSimplify

eq1FromOtherDemo=(m1+m2) L1 θ1''[t]+m2 L2 θ2''[t] Cos[θ1[t]-θ2[t]]+m2 L2 θ2'[t]^2 Sin[θ1[t]-θ2[t]]+g(m1+m2)Sin[θ1[t]]

(eq1-eq1FromOtherDemo)//FullSimplify

sol1=First@Solve[eq1==0,θ1''[t]]

Find equation of motion for θ2

eq2=D[D[L,θ2'[t]],t]-D[L,θ2[t]]+c θ2'[t]

eq2=eq2//FullSimplify

eq2FromOtherDemo=m2 L2 θ2''[t]+m2 L1 θ1''[t]Cos[θ1[t]-θ2[t]]-m2 L1 θ1'[t]^2 Sin[θ1[t]-θ2[t]]+m2 g Sin[θ2[t]]

(eq2-eq2FromOtherDemo)//FullSimplify

Expand[%]

sol2=First@Solve[eq2==0,θ2''[t]]

FullSimplify[eq2]

decouple the above 2 equations

sol=First@Solve[{eq1==0,eq2==0},{θ1''[t],θ2''[t]}]

Numerically solve the 2 equations found above for θ(t) and ϕ(t), use some arbitrary parameters

pars={L1->1,L2->1,m1->1,m2->1,g->9.8,c->1};
ic={θ1[0]==1.2  ,θ1'[0]==0,θ2[0]==3,θ2'[0]==0};
numericalSolution=First@NDSolve[Flatten[{eq1==0,eq2==0,ic}]/.pars,{θ1[t],θ2[t]},{t,0,10}]

frameLabel[t_]:={{None,None},{"t",t}}
at={Frame->True,ImageSize->250};

Row[{Plot[θ1[t]/.numericalSolution,{t,0,10},Evaluate@at,FrameLabel->frameLabel["θ1(t)"]],Plot[θ2[t]/.numericalSolution,{t,0,10},Evaluate@at,FrameLabel->frameLabel["θ2(t)"]]}]

Make plots of the solutions. Assume (Subscript[L, 1])=1,(Subscript[L, 2])=1,(Subscript[m, 1])=1,(Subscript[m, 2])=1,θ(0)=1,θ'(0)=0,ϕ(0)=3,ϕ'(0)=0

frameLabel[t_]:={{None,None},{"t",t}}
at={Frame->True,ImageSize->250};

Row[{Plot[θ1[t]/.numericalSolution,{t,0,10},Evaluate@at,FrameLabel->frameLabel["θ1(t)"]],Plot[θ2[t]/.numericalSolution,{t,0,10},Evaluate@at,FrameLabel->frameLabel["θ2(t)"]]}]

make phase portait plot

Plot[{θ1[t]/.numericalSolution,θ1'[t]/.numericalSolution},{t,0,10}]

Convert the above to state space form

sol=First@Solve[{eq1==0,eq2==0},{θ1''[t],θ2''[t]}];
secondOrderθ1=θ1''[t]/.sol
secondOrderθ2=θ2''[t]/.sol

repl={θ1[t]->x1[t],θ1'[t]->x2[t],θ2[t]->x3[t],θ2'[t]->x4[t]}

eq1=x1'[t]== x2[t]
eq2=x2'[t]==secondOrderθ1/.repl
eq3=x3'[t]== x4[t]
eq4=x4'[t]==secondOrderθ2/.repl

F=ma solution

Write the free body diagram

We have 4 unknowns (Subscript[T, 1]),(Subscript[T, 2]),θ(t),ϕ(t)

Write F=ma for forces radially along θ along for (Subscript[m, 1])

eq1=T1 -(Subscript[m, 1]) g Cos[θ[t]]-T2 Cos[ϕ[t]-θ[t]]== (Subscript[m, 1])((Subscript[L, 1]) θ'[t]^2);

Row[{"eq1 = ",eq1//FullSimplify//rep}]

Resolve forces perpendicular for (Subscript[m, 1])

eq2=T2 Sin[ϕ[t]-θ[t]] -(Subscript[m, 1]) g Sin[θ[t]]==(Subscript[m, 1])((Subscript[L, 1]) θ''[t]);

Row[{"eq2 = ",eq2//FullSimplify//rep}]

Write F=ma for forces radially along ϕ along for (Subscript[m, 2])

eq3=T2 -(Subscript[m, 2]) g Cos[ϕ[t]]==(Subscript[m, 2])((Subscript[L, 1]) θ'[t]^2 Cos[ϕ[t]-θ[t]]+(Subscript[L, 2]) ϕ'[t]^2-(Subscript[L, 1]) θ''[t]Sin[ϕ[t]-θ[t]]);

Row[{"eq3 = ",eq3//FullSimplify//rep}]

Resolve forces perpendicular for (Subscript[m, 2])

eq4=-(Subscript[m, 1]) g Sin[ϕ[t]]==(Subscript[m, 1])((Subscript[L, 2]) ϕ''[t]+(Subscript[L, 1]) θ''[t] Cos[ϕ[t]-θ[t]]+(Subscript[L, 1]) θ'[t]^2 Sin[ϕ[t]-θ[t]]);

Row[{"eq4 = ",eq4//FullSimplify//rep}]

Solve the above 4 equestions

solF=First@Solve[{eq1,eq2,eq3,eq4},{T1,T2,θ''[t],ϕ''[t]}];

Row[{"(Subscript[T, 1])=",T1/.solF//FullSimplify//rep}]
Row[{"(Subscript[T, 2])=",T2/.solF//FullSimplify//rep}]
Row[{"(Overscript[θ, ••])=",θ''[t]/.solF//FullSimplify//rep}]
Row[{"(Overscript[ϕ, ••])=",ϕ''[t]/.solF//FullSimplify//rep}]

We see that the solution from F=ma is the same as the solution obtained from Lagrange method.  To verify

(θ''[t]/.sol)-(θ''[t]/.solF)//Simplify

(ϕ''[t]/.sol)-(ϕ''[t]/.solF)//Simplify

We see from the above that the solutions are identical.

3 bob pendulum derivation

ke=(1/2)m1 (L1 θ1'[t])^2+(1/2)m2( (L1 θ1'[t] Cos[θ1[t]]+L2 θ2'[t] Cos[θ2[t]])^2+(L1 θ1'[t] Sin[θ1[t]]+L2 θ2'[t] Sin[θ2[t]])^2)+

(1/2)m3( (L1 θ1'[t] Cos[θ1[t]]+L2 θ2'[t] Cos[θ2[t]]+L3 θ3'[t] Cos[θ3[t]])^2+(L1 θ1'[t] Sin[θ1[t]]+L2 θ2'[t] Sin[θ2[t]]+L3 θ3'[t] Sin[θ3[t]])^2)

pe=m1 g( L1 (1-Cos[ θ1[t]]) )+m2 g(L1 (1-Cos[ θ1[t]]+ L2 (1-Cos[ θ2[t]])))+m3 g(L1 (1-Cos[ θ1[t]]+ L2 (1-Cos[ θ2[t]])+L3(1-Cos[θ3[t]])))

peOld=-L1 m1 g Cos[θ1[t]]-(L1 Cos[θ1[t]]+L2 Cos[θ2[t]])m2 g-(L1 Cos[θ1[t]]+L2 Cos[θ2[t]]+L3 Cos[θ3[t]])m3 g

lagOld=ke-peOld

lag=ke-pe

eq1=D[D[lag,θ1'[t]],t]-D[lag,θ1[t]]//Simplify

eq1Old=D[D[lagOld,θ1'[t]],t]-D[lagOld,θ1[t]]//Simplify

eq1-eq1Old//FullSimplify

eq2=D[D[lag,θ2'[t]],t]-D[lag,θ2[t]]

eq2=eq2/L2//FullSimplify

eq3=D[D[lag,θ3'[t]],t]-D[lag,θ3[t]]

eq3=eq3/L3//FullSimplify

try to decouple the 3 bob pendulum equations

sol=First@Solve[{eq1==0,eq2==0,eq3==0},{θ1''[t],θ2''[t],θ3''[t]}];
secondOrderθ1=θ1''[t]/.sol
secondOrderθ2=θ2''[t]/.sol
secondOrderθ3=θ3''[t]/.sol

repl={θ1[t]->x1[t],θ1'[t]->x2[t],θ2[t]->x3[t],θ2'[t]->x4[t],θ3[t]->x5[t],θ3'[t]->x6[t]}

eq1=x1'[t]== x2[t]
eq2=x2'[t]==secondOrderθ1/.repl
eq3=x3'[t]== x4[t]
eq4=x4'[t]==secondOrderθ2/.repl
eq5=x5'[t]== x6[t]
eq6=x6'[t]==secondOrderθ3/.repl

pars={L1->1,L2->1,L3->1,m1->1,m2->1,m3->1,g->9.8};
eq11=eq1/.pars
eq22=eq2/.pars
eq33=eq3/.pars
eq44=eq4/.pars
eq55=eq5/.pars
eq66=eq6/.pars

r=eq2/._==y_->y
Numerator[r]/.L1->0
Denominator[r]/.L1->0

ic={x1[0]==30 Degree  ,x2[0]==0,x3[0]==45 Degree  ,x4[0]==0,x5[0]==80 Degree  ,x6[0]==0};
sol=First@NDSolve[{eq1,eq2,eq3,eq4,eq5,eq6,x1[0]==30 Degree  ,x2[0]==0,x3[0]==45 Degree  ,x4[0]==0,x5[0]==80 Degree  ,x6[0]==0},{x1[t],x2[t],x3[t],x4[t],x5[t],x6[t]},{t,0,10}]

frameLabel[t_]:={{None,None},{"t",t}}
at={Frame->True,ImageSize->250};
Row[{Plot[x1[t]/.sol,{t,0,10},Evaluate@at,FrameLabel->frameLabel["x1(t)"]],Plot[x2[t]/.sol,{t,0,10},Evaluate@at,FrameLabel->frameLabel["x2(t)"]]}]

2 bob pendumul

ke=(1/2)m1 (L1 θ1'[t])^2+(1/2)m2( (L1 θ1'[t] Cos[θ1[t]]+L2 θ2'[t] Cos[θ2[t]])^2+(L1 θ1'[t] Sin[θ1[t]]+L2 θ2'[t] Sin[θ2[t]])^2)

pe=-L1 m1 g Cos[θ1[t]]-(L1 Cos[θ1[t]]+L2 Cos[θ2[t]])m2 g

lag=ke-pe

eq1=D[D[lag,θ1'[t]],t]-D[lag,θ1[t]]

eq11=g L1 m1 Sin[θ1[t]]+g L1 m2 Sin[θ1[t]]+g L1 m3 Sin[θ1[t]]-(1/2) m2 (-2 L1 Sin[θ1[t]] (θ1)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t])+2 L1 Cos[θ1[t]] (θ1)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t]))-(1/2) m3 (-2 L1 Sin[θ1[t]] (θ1)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t]+L3 Cos[θ3[t]] (θ3)'[t])+2 L1 Cos[θ1[t]] (θ1)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t]+L3 Sin[θ3[t]] (θ3)'[t]))+(L1^2) m1 (θ1)''[t]+(1/2) m2 (-2 L1 Sin[θ1[t]] (θ1)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t])+2 L1 Cos[θ1[t]] (θ1)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t])+2 L1 Cos[θ1[t]] (-L1 Sin[θ1[t]] (((θ1)'[t])^2)-L2 Sin[θ2[t]] (((θ2)'[t])^2)+L1 Cos[θ1[t]] (θ1)''[t]+L2 Cos[θ2[t]] (θ2)''[t])+2 L1 Sin[θ1[t]] (L1 Cos[θ1[t]] (((θ1)'[t])^2)+L2 Cos[θ2[t]] (((θ2)'[t])^2)+L1 Sin[θ1[t]] (θ1)''[t]+L2 Sin[θ2[t]] (θ2)''[t]))+(1/2) m3 (-2 L1 Sin[θ1[t]] (θ1)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t]+L3 Cos[θ3[t]] (θ3)'[t])+2 L1 Cos[θ1[t]] (θ1)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t]+L3 Sin[θ3[t]] (θ3)'[t])+2 L1 Cos[θ1[t]] (-L1 Sin[θ1[t]] (((θ1)'[t])^2)-L2 Sin[θ2[t]] (((θ2)'[t])^2)-L3 Sin[θ3[t]] (((θ3)'[t])^2)+L1 Cos[θ1[t]] (θ1)''[t]+L2 Cos[θ2[t]] (θ2)''[t]+L3 Cos[θ3[t]] (θ3)''[t])+2 L1 Sin[θ1[t]] (L1 Cos[θ1[t]] (((θ1)'[t])^2)+L2 Cos[θ2[t]] (((θ2)'[t])^2)+L3 Cos[θ3[t]] (((θ3)'[t])^2)+L1 Sin[θ1[t]] (θ1)''[t]+L2 Sin[θ2[t]] (θ2)''[t]+L3 Sin[θ3[t]] (θ3)''[t]))

eq11=eq11/.m3->0

eq11-eq1

eq2=D[D[lag,θ2'[t]],t]-D[lag,θ2[t]]

eq22=g L2 m2 Sin[θ2[t]]+g L2 m3 Sin[θ2[t]]-(1/2) m2 (-2 L2 Sin[θ2[t]] (θ2)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t])+2 L2 Cos[θ2[t]] (θ2)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t]))-(1/2) m3 (-2 L2 Sin[θ2[t]] (θ2)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t]+L3 Cos[θ3[t]] (θ3)'[t])+2 L2 Cos[θ2[t]] (θ2)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t]+L3 Sin[θ3[t]] (θ3)'[t]))+(1/2) m2 (-2 L2 Sin[θ2[t]] (θ2)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t])+2 L2 Cos[θ2[t]] (θ2)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t])+2 L2 Cos[θ2[t]] (-L1 Sin[θ1[t]] (((θ1)'[t])^2)-L2 Sin[θ2[t]] (((θ2)'[t])^2)+L1 Cos[θ1[t]] (θ1)''[t]+L2 Cos[θ2[t]] (θ2)''[t])+2 L2 Sin[θ2[t]] (L1 Cos[θ1[t]] (((θ1)'[t])^2)+L2 Cos[θ2[t]] (((θ2)'[t])^2)+L1 Sin[θ1[t]] (θ1)''[t]+L2 Sin[θ2[t]] (θ2)''[t]))+(1/2) m3 (-2 L2 Sin[θ2[t]] (θ2)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t]+L3 Cos[θ3[t]] (θ3)'[t])+2 L2 Cos[θ2[t]] (θ2)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t]+L3 Sin[θ3[t]] (θ3)'[t])+2 L2 Cos[θ2[t]] (-L1 Sin[θ1[t]] (((θ1)'[t])^2)-L2 Sin[θ2[t]] (((θ2)'[t])^2)-L3 Sin[θ3[t]] (((θ3)'[t])^2)+L1 Cos[θ1[t]] (θ1)''[t]+L2 Cos[θ2[t]] (θ2)''[t]+L3 Cos[θ3[t]] (θ3)''[t])+2 L2 Sin[θ2[t]] (L1 Cos[θ1[t]] (((θ1)'[t])^2)+L2 Cos[θ2[t]] (((θ2)'[t])^2)+L3 Cos[θ3[t]] (((θ3)'[t])^2)+L1 Sin[θ1[t]] (θ1)''[t]+L2 Sin[θ2[t]] (θ2)''[t]+L3 Sin[θ3[t]] (θ3)''[t]));
eq22=eq22/.m3->0

eq2-eq22

pars={L1->1,L2->1,m1->1,m2->1,g->9.8};
ic={θ1[0]==1.2  ,θ1'[0]==0,θ2[0]==1.2  ,θ2'[0]==0};
numericalSolution=First@NDSolve[Flatten[{eq1==0,eq2==0,ic}]/.pars,{θ1[t],θ2[t]},{t,0,10}]

NDSolve[%,{θ1[t],θ2[t]},{t,0,10}]

eq33=g L3 m3 Sin[θ3[t]]-(1/2) m3 (-2 L3 Sin[θ3[t]] (θ3)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t]+L3 Cos[θ3[t]] (θ3)'[t])+2 L3 Cos[θ3[t]] (θ3)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t]+L3 Sin[θ3[t]] (θ3)'[t]))+(1/2) m3 (-2 L3 Sin[θ3[t]] (θ3)'[t] (L1 Cos[θ1[t]] (θ1)'[t]+L2 Cos[θ2[t]] (θ2)'[t]+L3 Cos[θ3[t]] (θ3)'[t])+2 L3 Cos[θ3[t]] (θ3)'[t] (L1 Sin[θ1[t]] (θ1)'[t]+L2 Sin[θ2[t]] (θ2)'[t]+L3 Sin[θ3[t]] (θ3)'[t])+2 L3 Cos[θ3[t]] (-L1 Sin[θ1[t]] (((θ1)'[t])^2)-L2 Sin[θ2[t]] (((θ2)'[t])^2)-L3 Sin[θ3[t]] (((θ3)'[t])^2)+L1 Cos[θ1[t]] (θ1)''[t]+L2 Cos[θ2[t]] (θ2)''[t]+L3 Cos[θ3[t]] (θ3)''[t])+2 L3 Sin[θ3[t]] (L1 Cos[θ1[t]] (((θ1)'[t])^2)+L2 Cos[θ2[t]] (((θ2)'[t])^2)+L3 Cos[θ3[t]] (((θ3)'[t])^2)+L1 Sin[θ1[t]] (θ1)''[t]+L2 Sin[θ2[t]] (θ2)''[t]+L3 Sin[θ3[t]] (θ3)''[t]));
eq33/.m3->0

3 bobs pendulum equations

ode1

ode2

ode3