Application Center - Maplesoft

# Dynamics of a double pendulum

You can switch back to the summary page by clicking here.

doublependulum.mws

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
```

```Warning, the name changecoords has been redefined
```

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:

> od:

> for i from 1 by 1 to 4 do

> f0||i:=f0[num[i]];

> DL||i:=DL[num[i]];

> o||i:=o[num[i]];

> od:

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

> for i from 1 by 1 to 4 do

> for j from 1 by 1 to 4 do

> T[j,i]:=ev[i][3][1][j];

> od;

> ewk[i]:=ev[i][1]:

> od:

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:

> if (i<>j) then AS[i,j]:=0 fi:

> od:

> od:

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:

Displacements

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

Velocities

> 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:

> od:

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

> od:

> 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.