Dynamics of a Double Pendulum
2000 Harald Kammerer, GERB Schwingungsisolierungen, Germany http://www.gerb.com
This worksheet shows the calculation of the motion of a double pendulum
The author expects that this worksheet will only be used for teaching and educational purposes and not for commerical profit without contacting the author for a licensed agreement.
> restart;
Introduction
Consider a double pendulum system. At first the eigenfrequencies and the eigenmodes of the system are calculated. Second the motion of the system caused by specific initial coditions is calculated.
Index 1 belongs to the upper mass (m1, l1, d1, w1) , index 2 to the lower mass (m2, l2, d2, w2).
All units are kg, m, N and sec
The calculation is valid for small displacements, so that sin(x)=x and cos(x)=1.
Initialisation
> with(linalg):with(plots):with(plottools):with(linalg):
Warning, the name changecoords has been redefined
Warning, the protected names norm and trace have been redefined and unprotected
Description of the Mechanical Modell
Geometry
upper mass
> m1:=10.:l1:=1.:d1:=1.:
lower mass
> m2:=10.:l2:=1.:d2:=1.5:
weight force
> G1:=m1*9.81:
> G2:=m2*9.81:
Initial Conditions
> w10:=0:
> w20:=0:
> wp10:=0:
> wp20:=1:
System Matrices
Substitution of the derivations
>
System of Equations of Motion
Consider the momentum around point A.
Note: G*sin(phi)=G*phi; G*cos(phi)=G.
Mass-, Damping- Stiffness- and System Matrices
Mass
> MM:=matrix(2,2,[m1,0,0,m2]);
Damping
> DM:=matrix(2,2,[d1,0,0,d2]);
Stiffness
Collect the coefficients of the displacement terms of the system equation
> k11:=-coeff(rhs(eq1),w1(t)):
> k12:=-coeff(rhs(eq1),w2(t)):
> k21:=-coeff(rhs(eq2),w1(t)):
> k22:=-coeff(rhs(eq2),w2(t)):
> KM:=matrix(2,2,[k11,k12,k21,k22]);
System Matrix
Transform the 2 homogen equation of motion of degree 2 into 4 differencial equation of degree 1
> smd11:=matrix(2,2,[0,0,0,0]):smd12:=matrix(2,2,[1,0,0,1]):
> smd21:=-multiply(inverse(MM),KM):smd22:=-multiply(inverse(MM),DM):
> SMD:=matrix(4,4):
> SMD:=copyinto(smd11,SMD,1,1):SMD:=copyinto(smd12,SMD,1,3):SMD:=copyinto(smd21,SMD,3,1):SMD:=copyinto(smd22,SMD,3,3);
Eigenvalues and Eigenmodes
Solve the Eigenvalue Problem
> ewd:=eigenvals(SMD);
Eigenfrequency and Damping of the Eigenmodes
Eigenfrequency: imaginary part of the eigenvalue devided by 2 Pi.
Damping of the eigenmode: real part of the eigenvalue devided by the amount.
Note: only valid if the damping is undercritical.
> f01:=evalf(Im(ewd[1])/2/Pi):f02:=evalf(Im(ewd[2])/2/Pi):f03:=evalf(Im(ewd[3])/2/Pi):f04:=evalf(Im(ewd[4])/2/Pi):
> o1:=sqrt(Im(ewd[1])**2+Re(ewd[1])**2):o2:=sqrt(Im(ewd[2])**2+Re(ewd[2])**2):o3:=sqrt(Im(ewd[3])**2+Re(ewd[3])**2):o4:=sqrt(Im(ewd[4])**2+Re(ewd[4])**2):
> D1:=-Re(ewd[1])/o1:D2:=-Re(ewd[2])/o2:D3:=-Re(ewd[3])/o3:D4:=-Re(ewd[4])/o4:
Sort according the Eigenvalues
> o:=[o1,o2,o3,o4]:f0:=[f01,f02,f03,f04]:DL:=[D1,D2,D3,D4]:
> lstsrt:=sort([o[1],o[2],o[3],o[4]]):
> for i from 1 by 1 to 4 do
> for j from 1 by 1 to 4 do
> if (lstsrt[i]=o[j]) then
> num[i]:=j:
> if (i>1 and num[i]<>num[i-1]) then j:=4 fi:
> fi:
> od:
> f0||i:=f0[num[i]];
> DL||i:=DL[num[i]];
> o||i:=o[num[i]];
View the Eigenvalues Numerical
(Eigenvalues with positiv and negativ sign but the same amount belong together)
> f01;DL1;
> f02;DL2;
> f03;DL3;
> f04;DL4;
Eigenmodes, Transformation Matrix
The eigenvectors arranged in form of a matrix yields the transformation matrix T to calculate the Motion in main coordinates
> ev:=eigenvects(SMD):
> unassign('i','j');
> T:=matrix(4,4):
> T[j,i]:=ev[i][3][1][j];
> od;
> ewk[i]:=ev[i][1]:
Transform System Matrixs into Main Coordinates
> AS:=multiply(inverse(T),multiply(SMD,T)):
elmininate numerical inaccuraties
> for i from 1 by 1 to 4 do for j from 1 by 1 to 4 do
> if (i<>j) then AS[i,j]:=0 fi:
Calculation of the Motion
Transform the Initial Conditions into Main Coordinates
Vector of initial conditions
> w0:=vector(4,[w10,w20,wp10,wp20]);
Transformation in main coordinates
> z0:=multiply(inverse(T),w0);
> AB1:=z1(0)=z0[1];AB2:=z2(0)=z0[2];AB3:=z3(0)=z0[3];AB4:=z4(0)=z0[4];
Main coordinates and derivations
> Z:=matrix(4,1,[z1(t),z2(t),z3(t),z4(t)]);
> Zp:=matrix(4,1,[diff(z1(t),t),diff(z2(t),t),diff(z3(t),t),diff(z4(t),t)]);
Equation of Motion in Main Coordinates
> BGL:=Zp=AS&*Z;
separated equations
> bgl1:=lhs(evalm(BGL))[1,1]=rhs(evalm(BGL))[1,1];
> bgl2:=lhs(evalm(BGL))[2,1]=rhs(evalm(BGL))[2,1];
> bgl3:=lhs(evalm(BGL))[3,1]=rhs(evalm(BGL))[3,1];
> bgl4:=lhs(evalm(BGL))[4,1]=rhs(evalm(BGL))[4,1];
Solution
> lsg:=dsolve({bgl1,bgl2,bgl3,bgl4,AB1,AB2,AB3,AB4},{z1(t),z2(t),z3(t),z4(t)}):
> assign(lsg);
Transformation in Real Coordinates
> ZL:=matrix(4,1,[combine(expand(evalf(z1(t)))),combine(expand(evalf(z2(t)))),combine(expand(evalf(z3(t)))),combine(expand(evalf(z4(t))))]):
> YL:=multiply(T,ZL):
Displacements
> w1(t):=simplify(expand(YL[1,1])):
> w2(t):=simplify(expand(YL[2,1])):
Velocities
> v1(t):=simplify(expand(YL[3,1])):
> v2(t):=simplify(expand(YL[4,1])):
Graphical View of the Time History
Number of time steps
> ndt:=100:
Timestep
> dt:=0.5:
> P1:=plot(Re(w1(t)),t=0..ndt*dt,color=red):
> P2:=plot(Re(w2(t)),t=0..ndt*dt,color=blue):
> display({P1,P2});
> PV1:=plot(Re(v1(t)),t=0..ndt*dt,color=red):
> PV2:=plot(Re(v2(t)),t=0..ndt*dt,color=blue):
> display({PV1,PV2});
Phase Curve
> PPh1:=plot([Re(w1(t)),Re(v1(t)),t=0..ndt*dt],color=red):
> PPh2:=plot([Re(w2(t)),Re(v2(t)),t=0..ndt*dt],color=blue):
> display({PPh1,PPh2});
Animation
> FIX1:=curve([[-0.2,0],[-0.05,0]],color=black,thickness=2):
> FIX2:=circle([0,0],0.05,color=black):
> FIX3:=curve([[0.05,0],[0.2,0]],color=black,thickness=2):
> FIX:=display([FIX1,FIX2,FIX3]):
> MAN01:=circle([-0.2,-(l1+l2)+0.05],0.05,color=black,thickness=2):
> MAN02:=curve([[-0.2,-(l1+l2)],[-0.2,-(l1+l2)-0.15]],color=black,thickness=2):
> MAN03:=curve([[-0.2,-(l1+l2)-0.05],[-0.25,-(l1+l2)]],color=black,thickness=2):
> MAN04:=curve([[-0.2,-(l1+l2)-0.05],[-0.15,-(l1+l2)]],color=black,thickness=2):
> MAN05:=curve([[-0.2,-(l1+l2)-0.15],[-0.25,-(l1+l2)-0.3]],color=black,thickness=2):
> MAN06:=curve([[-0.2,-(l1+l2)-0.15],[-0.15,-(l1+l2)-0.3]],color=black,thickness=2):
> MAN0:=display([MAN01,MAN02,MAN03,MAN04,MAN05,MAN06]):
> MAN11:=circle([-0.2,-(l1+l2)+0.05],0.05,color=black):
> MAN12:=curve([[-0.2,-(l1+l2)],[-0.2,-(l1+l2)-0.15]],color=black,thickness=2):
> MAN13:=curve([[-0.2,-(l1+l2)-0.05],[-0.25,-(l1+l2)]],color=black,thickness=2):
> MAN14:=curve([[-0.2,-(l1+l2)-0.05],[-0.1,-(l1+l2)]],color=black,thickness=2):
> MAN15:=curve([[-0.2,-(l1+l2)-0.15],[-0.25,-(l1+l2)-0.3]],color=black,thickness=2):
> MAN16:=curve([[-0.2,-(l1+l2)-0.15],[-0.15,-(l1+l2)-0.3]],color=black,thickness=2):
> MAN1:=display([MAN11,MAN12,MAN13,MAN14,MAN15,MAN16]):
> MAN21:=circle([-0.3,-(l1+l2)],0.05,color=black):
> MAN22:=curve([[-0.3,-(l1+l2)-0.05],[-0.3,-(l1+l2)-0.20]],color=black,thickness=2):
> MAN23:=curve([[-0.3,-(l1+l2)-0.10],[-0.35,-(l1+l2)-0.05]],color=black,thickness=2):
> MAN24:=curve([[-0.3,-(l1+l2)-0.10],[-0.25,-(l1+l2)-0.05]],color=black,thickness=2):
> MAN25:=curve([[-0.3,-(l1+l2)-0.20],[-0.25,-(l1+l2)-0.15]],color=black,thickness=2):
> MAN26:=curve([[-0.3,-(l1+l2)-0.20],[-0.25,-(l1+l2)-0.18]],color=black,thickness=2):
> MAN2:=display([MAN21,MAN22,MAN23,MAN24,MAN25,MAN26]):
> MAN[0]:=MAN0:
> MAN[1]:=MAN1:
> MAN[2]:=MAN1:
> for i from 3 by 1 to ndt-1 do
> MAN[i]:=MAN2:
> MAN[ndt]:=MAN0:
> for i from 0 by 1 to ndt do
> x1:=evalf(Re(subs(t=i*dt,w1(t)))):
> x2:=evalf(Re(subs(t=i*dt,w2(t)))):
> PLI[i]:=curve([[0,0],[x1,-l1],[x2,-(l1+l2)]]):
> MASS1[i]:=disk([x1,-l1], 0.1, color=red):
> MASS2[i]:=disk([x2,-(l1+l2)], 0.1, color=blue):
> ANIM[i]:=display({PLI[i],MASS1[i],MASS2[i],FIX,MAN[i]});
> display([seq(ANIM[i],i=0..ndt)],insequence=true,scaling=constrained,axes=none,title=`Double Pendulum`);
Disclaimer: While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.
