Modal Analysis for 3 pendulum with springs problem

by Nasser M Abbasi

This small note computes the eigenvectors of a 3 degrees freedom by modal analysis and without modal analysis and analysis the difference between the aproaches.

This digram below describes the problem. We use Lagrangian formulation to determine the equation of motions, then use modal analysis to decouple the system and solve it. In this system , the springs are attached at a distant α From the edge. Each pendulum has length L and has masses m1,m2,m3 attached to the end. To obtain a numerical solution, we assume some initial conditions such as  θ(0)={Pi/4,Pi/4,Pi/4} and (OverDot[θ])={0,0,0}

Modal analysis by decoupling

First define some notations to use

Needs["Notation`"]

Symbolize[(NotationTemplateTag[(Subscript[θ, 1])])];
Symbolize[(NotationTemplateTag[(Subscript[θ, 2])])];
Symbolize[(NotationTemplateTag[(Subscript[θ, 3])])];
Symbolize[(NotationTemplateTag[(Subscript[r, 1])])];
Symbolize[(NotationTemplateTag[(Subscript[r, 2])])];
Symbolize[(NotationTemplateTag[(Subscript[r, 3])])];
Symbolize[(NotationTemplateTag[(OverTilde[K])])];

Define a function which accepts the kinetic and potential energy and return back the stiffness and the mass matrix

obtainStiffnessAndMassMatrix[ke_,u_]:=Module[{lagrangian,eq1,eq2,eq3,meq1,meq2,meq3,seq1,seq2,seq3,stiffnessMatrix,massMatrix,any},
lagrangian=ke-u;
eq1=D[D[lagrangian,(Subscript[θ, 1])'[t]],t]-D[lagrangian,(Subscript[θ, 1])[t]]==0;
eq2=D[D[lagrangian,(Subscript[θ, 2])'[t]],t]-D[lagrangian,(Subscript[θ, 2])[t]]==0;
eq3=D[D[lagrangian,(Subscript[θ, 3])'[t]],t]-D[lagrangian,(Subscript[θ, 3])[t]]==0;

(*linearize using small angle approximation*)
eq1=eq1/.Sin[(Subscript[θ, 1])[t]]->(Subscript[θ, 1])[t];
eq2=eq2/.Sin[(Subscript[θ, 2])[t]]->(Subscript[θ, 2])[t];
eq3=eq3/.Sin[(Subscript[θ, 3])[t]]->(Subscript[θ, 3])[t];

Print["Equation of motion are"];
Print[{eq1,eq2,eq3}//TableForm];

seq1=eq1/. (Subscript[θ, 1])''[t]_ ->0;
seq2=eq2/. (Subscript[θ, 2])''[t]_ ->0;
seq3=eq3/. (Subscript[θ, 3])''[t]_ ->0;

z1=eq1/.any_==0->any;
z2=seq1/.any_==0->any;
meq1= (eq1/.any_==0->any)-(seq1/.any_==0->any)//Simplify;

z1=eq2/.any_==0->any;
z2=seq2/.any_==0->any;
meq2=(z1-z2)//FullSimplify;

meq3=(eq3/.any_==0->any)-(seq3/.any_==0->any)//Simplify;

r=Normal[CoefficientArrays[{seq1,seq2,seq3},{(Subscript[θ, 1])[t],(Subscript[θ, 2])[t],(Subscript[θ, 3])[t]}]];
stiffnessMatrix=Collect[r[[2,All]],α];
r=Normal[CoefficientArrays[{meq1,meq2,meq3},{(Subscript[θ, 1])''[t],(Subscript[θ, 2])''[t],(Subscript[θ, 3])''[t]}]];
massMatrix=Collect[r[[2,All]],α];
{massMatrix,stiffnessMatrix}
]

Now define the kinetic and potenatial energy

ke=(1/2)(Subscript[m, 1]) (((L (Subscript[θ, 1])'[t]))^2)+(1/2)(Subscript[m, 2]) (((L (Subscript[θ, 2])'[t]))^2)+(1/2)(Subscript[m, 3]) (((L (Subscript[θ, 3])'[t]))^2);
u=(1/2) (Subscript[k, 1]) ((( α (Subscript[θ, 2])[t]-α (Subscript[θ, 1])[t]))^2)+(1/2) (Subscript[k, 2]) ((( α (Subscript[θ, 3])[t]-α (Subscript[θ, 2])[t]))^2)+(Subscript[m, 1]) g  L(1-Cos[ (Subscript[θ, 1])[t]])+(Subscript[m, 2]) g  L(1-Cos[ (Subscript[θ, 2])[t]])+(Subscript[m, 3]) g L(1-Cos[ (Subscript[θ, 3])[t]]);

Now call the above function to generate the stiffness and mass matrix. It also prints the 3 equations of motion

{massMatrix,stiffnessMatrix}=obtainStiffnessAndMassMatrix[ke,u];

Now print the STIFFNESS and MASS matrix

vars={(Subscript[θ, 1])''[t],(Subscript[θ, 2])''[t],(Subscript[θ, 3])''[t]};
deps={(Subscript[θ, 1])[t],(Subscript[θ, 2])[t],(Subscript[θ, 3])[t]};
Print[MatrixForm[massMatrix],MatrixForm[Transpose[List[vars]]],"+",
MatrixForm[stiffnessMatrix],MatrixForm[Transpose[List[deps]]]];

Now that we have the stiffness and mass matrix, we can perform modal analysis. Start by doing the first transformation

invMassMatrix=MatrixPower[massMatrix,((-1)/2)];
Print["(OverTilde[K])= ",MatrixForm[(OverTilde[K])=invMassMatrix .stiffnessMatrix .invMassMatrix]];

Now define some numerical values to use for the rest of the analysis amd generate (OverTilde[K])

values={(Subscript[m, 1])->1,(Subscript[m, 2])->2,(Subscript[m, 3])->3,(Subscript[k, 1])->10,L->1,α->.5,(Subscript[k, 2])->2,g->9.8};
invMassMatrix=invMassMatrix/.values;
massMatrix=massMatrix/.values;
MatrixForm[(OverTilde[K])=(OverTilde[K])/.values]

Find  eigenvalues of the (OverTilde[K]) matrix

{lambda,eigvect}=Eigensystem[(OverTilde[K])];
MatrixForm[lambda]

Find  eigenvectors  of the (OverTilde[K]) matrix

eigvect=Transpose[eigvect];
MatrixForm[eigvect=eigvect/Map[Norm,eigvect]]

Now find Λ matrix

MatrixForm[Λ=Transpose[eigvect].(OverTilde[K]).eigvect//Chop]

Hence the decoupled system of differential equations is

vars={(Subscript[r, 1])''[t],(Subscript[r, 2])''[t],(Subscript[r, 3])''[t]};
deps={(Subscript[r, 1])[t],(Subscript[r, 2])[t],(Subscript[r, 3])[t]};
Print[MatrixForm[Transpose[List[vars]]],"+",
MatrixForm[Λ],MatrixForm[Transpose[List[deps]]]];

Now convert the IC from θ (t) space to r (t) space

initialθ={Pi/4,Pi/2,Pi/4};initialR=Transpose[eigvect].MatrixPower[massMatrix,(1/2)].initialθ//Chop

initialθDot={0,0,0};
initialRDot=Transpose[eigvect].MatrixPower[massMatrix,(1/2)].initialθDot

Now solve the r (t) system

eq=(Subscript[r, 1])''[t]+Λ[[1,1]](Subscript[r, 1])[t]==0;
solr1=First[DSolve[{eq,(Subscript[r, 1])[0]==initialR[[1]],(Subscript[r, 1])'[0]==initialRDot[[1]]},(Subscript[r, 1])[t],t]];
Print["(Subscript[r, 1])[t]=",solr1=(Subscript[r, 1])[t]/.solr1]

eq=(Subscript[r, 2])''[t]+Λ[[2,2]](Subscript[r, 2])[t]==0;
solr2=First[DSolve[{eq,(Subscript[r, 2])[0]==initialR[[2]],(Subscript[r, 2])'[0]==initialRDot[[2]]},(Subscript[r, 2])[t],t]];
Print["(Subscript[r, 2])[t]=",solr2=(Subscript[r, 2])[t]/.solr2]

eq=(Subscript[r, 3])''[t]+Λ[[3,3]](Subscript[r, 3])[t]==0;
solr3=First[DSolve[{eq,(Subscript[r, 3])[0]==initialR[[3]],(Subscript[r, 3])'[0]==initialRDot[[3]]},(Subscript[r, 3])[t],t]];
Print["(Subscript[r, 3])[t]=",solr3=(Subscript[r, 3])[t]/.solr3]

Now convert solution from r (t) to θ (t)

MatrixForm[solθ=invMassMatrix.eigvect.{{solr1},{solr2},{solr3}}]

solθ1=solθ[[1]];
solθ2=solθ[[2]];
solθ3=solθ[[3]];

Now plot the solutions

Print["(Subscript[θ, 1])[t]=",solθ1=First[Simplify[Chop[ExpToTrig[solθ1/.values]]]]];

Print["(Subscript[θ, 2])[t]=",solθ2=First[FullSimplify[Chop[ExpToTrig[solθ2/.values]]]]];

Print["(Subscript[θ, 3])[t]=",solθ3=First[Simplify[Chop[ExpToTrig[solθ3/.values]]]]];

Plot[{solθ1,solθ2,solθ3},{t,0,10}]

Modal analysis without decoupling

massMatrix

stiffnessMatrix