DetailedDualDoublePenDerivation.html

This worksheet demonstrates the Euler Lagrange method of deriving the equations of motion of a dual double pendulum system, that is, a pendulum with one upper arm and two lower arms. We seek to extend this construction to n lower arms. We assume massless arms with mass concentrated at the ends, which we call the "bobs"

First, we give the x and y coordinate positions of the bobs.

> Xone:=l1*sin(theta1(t));

> Yone:=l1*cos(theta1(t));

> Xtwo:=l2*sin(theta2(t))+Xone;

> Ytwo:=l2*cos(theta2(t))+Yone;

> Xthree:=l3*sin(theta3(t))+Xone;

> Ythree:=l3*cos(theta3(t))+Yone;

Now we take the derivatives of these with respect to t:

> X1dot:=diff(Xone,t);

> Y1dot:=diff(Yone,t);

> X2dot:=diff(Xtwo,t);

> Y2dot:=diff(Ytwo,t);

> X3dot:=diff(Xthree,t);

> Y3dot:=diff(Ythree,t);

Here we find the kinetic energy:

> Tone:= 1/2*m1*(X1dot^2 + Y1dot^2);

> Ttwo:= 1/2*m2*(X2dot^2 + Y2dot^2);

> Tthree:= 1/2*m3*(X3dot^2 + Y3dot^2);

And the potential energy:

> Uone:= -1*m1*g*l1*cos(theta1(t));

> Utwo:= -1*m2*g*(cos(theta1(t))*l1 + cos(theta2(t))*l2);

> Uthree:= -1*m3*g*(cos(theta1(t))*l1 + cos(theta3(t))*l3);

We form the Lagrangian:

> L:= simplify((Tone+Ttwo+Tthree)-(Uone+Utwo+Uthree));

In order to use the Euler-Lagrange equation

we must take derivatives with respect to compound variables like diff(theta1(t),t), and Maple cannot take these derivatives. We must substitute in variables for these expressions, take the derivatives then substitute back.

> dLdtheta1:=simplify(subs(u1=theta1(t),diff(subs(theta1(t)=u1,L),u1)));
dLdtheta1dot:=simplify(subs(v1=diff(theta1(t),t),diff(subs(diff(theta1(t),t)=v1,L),v1)));

> dLdtheta2:=simplify(subs(u2=theta2(t),diff(subs(theta2(t)=u2,L),u2)));
dLdtheta2dot:=simplify(subs(v2=diff(theta2(t),t),diff(subs(diff(theta2(t),t)=v2,L),v2)));

> dLdtheta3:=simplify(subs(u3=theta3(t),diff(subs(theta3(t)=u3,L),u3)));
dLdtheta3dot:=simplify(subs(v3=diff(theta3(t),t),diff(subs(diff(theta3(t),t)=v3,L),v3)));

Here we take the time derivative of the derivative of the Lagrangian with respect to theta_i dot

> dotdLdtheta1dot:=simplify(diff(dLdtheta1dot,t));
dotdLdtheta2dot:=simplify(diff(dLdtheta2dot,t));
dotdLdtheta3dot:=simplify(diff(dLdtheta3dot,t));

The Euler-Lagrange equations are assembled:

> ELtheta1:=simplify(dLdtheta1-dotdLdtheta1dot);
ELtheta2:=simplify(dLdtheta2-dotdLdtheta2dot);
ELtheta3:=simplify(dLdtheta3-dotdLdtheta3dot);

Here for the sake of easy simulation, we set most of the variables to one:

> ELtheta1a:=subs({g=9.8,m1=1,m2=1,m3=1,l1=1,l2=1,l3=1},ELtheta1);
ELtheta2a:=subs({g=9.8,m1=1,m2=1,m3=1,l1=1,l2=1,l3=1},ELtheta2);
ELtheta3a:=subs({g=9.8,m1=1,m2=1,m3=1,l1=1,l2=1,l3=1},ELtheta3);

Here we put together some initial conditions:

> with(DEtools):epsilon:=10^(-1);IC:=theta1(0)=Pi,D(theta1)(0)=-epsilon,theta2(0)=Pi,D(theta2)(0)=0,theta3(0)=Pi,D(theta3)(0)=epsilon;

The function 'sol' is the simultaneous solution to the three Euler-Lagrange equations.

> sol:=dsolve({ELtheta1a,ELtheta2a,ELtheta3a,IC},{theta1(t),theta2(t),theta3(t)},numeric,method=lsode,output=listprocedure):

Now we give the names θ1a(t), θ2a(t), θ3a(t), for the variables θ1(t), θ2(t), θ3(t), with the numerical solution for these variables substituted for θ1(t), θ2(t), θ3(t):

> theta1a:=subs(sol,theta1(t)): theta2a:=subs(sol,theta2(t)): theta3a:=subs(sol,theta3(t)):

And we plot a short time series:

> with(plots):plot([theta1a(t),theta2a(t),theta3a(t)],t=0..3);

These functions x1, x2, x3, y1, y2, y3 are the positions of the ends of the bobs:

> x1:=t->sin(theta1a(t)):

> y1:=t->-cos(theta1a(t)):

> x2:=t->x1(t)+sin(theta2a(t)):

> y2:=t->y1(t)-cos(theta2a(t)):

> x3:=t->x1(t)+sin(theta3a(t)):

> y3:=t->y1(t)-cos(theta3a(t)):

The animation will consist of lines whose endpoints correspond to the positions of the bobs:

> PL:=t->plot([[0,0],[x1(t),y1(t)],[x2(t),y2(t)],[x1(t),y1(t)],[x3(t),y3(t)]],style=line,thickness=2,color=black,title=cat("time= ",convert(evalf(t,3),string))):

> PLa:=seq(PL(i/10),i=1..300):

> display(PLa,insequence=true,scaling=constrained,axes=boxed);

> ELt1:=subs({theta1(t)=u1,theta2(t)=u2,theta3(t)=u3,diff(theta1(t),t)=v1,diff(theta2(t),t)=v2,diff(theta3(t),t)=v3,diff(theta1(t),t,t)=w1,diff(theta2(t),t,t)=w2,diff(theta3(t),t,t)=w3},ELtheta1);
ELt2:=subs({theta1(t)=u1,theta2(t)=u2,theta3(t)=u3,diff(theta1(t),t)=v1,diff(theta2(t),t)=v2,diff(theta3(t),t)=v3,diff(theta1(t),t,t)=w1,diff(theta2(t),t,t)=w2,diff(theta3(t),t,t)=w3},ELtheta2);
ELt3:=subs({theta1(t)=u1,theta2(t)=u2,theta3(t)=u3,diff(theta1(t),t)=v1,diff(theta2(t),t)=v2,diff(theta3(t),t)=v3,diff(theta1(t),t,t)=w1,diff(theta2(t),t,t)=w2,diff(theta3(t),t,t)=w3},ELtheta3);