Problem #1)

Clear["Global`*"]

T=(1/2)m b^2 (dθ1^2+dθ2^2+dθ3^2)

V=(1/2)k b^2 ((θ1-θ2)^2+(θ2-θ3)^2+(θ3-θ1)^2) //Expand

lag=T-V

D[T,dθ1]

Tmat=m b^2 DiagonalMatrix[{1,1,1}];
Tmat //MatrixForm

D[V,θ1]

Vmat=k b^2{{2,-1,-1},{-1,2,-1},{-1,-1,2}};
Vmat //MatrixForm

mat= Vmat-Tmat ω2;
mat //MatrixForm

d[D[lag,dθ1],t] - D[lag,θ1]

d[D[lag,dθ2],t] - D[lag,θ2]

d[D[lag,dθ3],t] - D[lag,θ3]

sol=Solve[Det[mat] ==0,ω2]

eq1=mat .{a,b,c} ==0 //Thread;
eq1 //Column

norm = a^2+b^2+c^2==1

eq2=Join[eq1,{norm}];
eq2 //Column

sol1=Solve[eq2 /.sol[[1]],{a,b,c}] //Last

ev1={a,b,c} /.sol1

sol2=Solve[eq2 /.sol[[2]],{a,b,c}]

(* Let's look for a solution where one of them is at rest *)
ev2={a,b,c} /.Last[sol2] /.a->0

ev3=Cross[ev1,ev2]

eVecs={ev1,ev2,ev3}

Part B Look at motion

values={k->1,m->1};

sol

mode1=ev1 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[1]]

mode2=ev2 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[2]]

mode3=ev3 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[3]]

Plot[mode1 +{-2,0,2}/.values //Re,{t,0,10}]

Plot[mode2 +{-2,0,2}/.values //Re,{t,0,10}]

Plot[mode3 +{-2,0,2}/.values //Re,{t,0,10}]

Plot[mode2 + mode3 +{-2,0,2}/.values //Re,{t,0,10}]

(* a cross check  *)
Transpose[eVecs]==Inverse[eVecs]
eVecs.Transpose[eVecs]

Vdiag=eVecs.Vmat.Transpose[eVecs] //Simplify ;
Vdiag//MatrixForm

Tdiag=eVecs.Tmat.Transpose[eVecs] //Simplify;
Tdiag //MatrixForm

Problem #2)

Clear["Global`*"]

T=(1/2)m v2+(1/2) Ibar ωBar^2;
ωBar=(((x1'[t]-x2'[t]))/L);
Ibar=(1/12) m L^2;
vel=(((x1'[t]+x2'[t]))/2);
v2=(vel^2)
T//Expand//Factor

V=-m g (((x1[t]+x2[t]))/2)+(1/2)k (x1[t]^2+x2[t]^2)//Simplify

lag=T-V;

D[D[T,x1'[t]],t] //Expand

D[D[T,x2'[t]],t] //Expand

eq1= D[D[lag,x1'[t]],t] - D[lag,x1[t]] //Expand

eq2= D[D[lag,x2'[t]],t] - D[lag,x2[t]]//Expand

Tmat=(m/6){{2,1},{1,2}};
Tmat //MatrixForm

D[V,x1[t]] //Expand

D[V,x2[t]] //Expand

Vmat=k{{1,0},{0,1}};

Vmat //MatrixForm

Tmat //MatrixForm

mat= Vmat-Tmat ω2;
mat//MatrixForm

sol=Solve[Det[mat] ==0,ω2]

eq0=mat .{a,b} ==0 //Thread

norm = a^2+b^2==1

eqs=Join[eq0,{norm}]

sol1=Solve[eqs /.sol[[1]],{a,b}]

ev1={a,b} /.sol1[[1]]

sol2=Solve[eqs /.sol[[2]],{a,b}]

ev2={a,b} /.sol2[[1]]

ev1.ev2

eVecs={ev1,ev2}

eVecs.Vmat.Transpose[eVecs] //Simplify //MatrixForm

eVecs.Tmat.Transpose[eVecs] //Simplify //MatrixForm

Part B Look at motion

values={k->1,m->1};

sol

mode1=ev1 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[1]]

mode2=ev2 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[1]]

mode1 +{-1,1} /.values

Plot[mode1 +{-1,1}/.values //Re,{t,0,10}]

Plot[mode2 +{-1,1}/.values //Re,{t,0,10}]

Plot[(1/3)mode1 + (2/3)mode2  +{-2,2}/.values //Re,{t,0,10}]

{x1,x2}==eVecs.{n1,n2}

{x1,x2}==eVecs.{n1,n2} /.{n1->1/2,n2->1/2}

mix[a_]:=a mode1 +(1-a)mode2 /.values //Re
mix[1/2] /.{t->0} //Simplify

doPlot[a_:1/2]:=Plot[a (mode1 +{-1,1}/.values)+(1-a)(mode2 +{-1,1}/.values)
 //Re,{t,0,10}]

doPlot[1/4]

Problem #3)

Clear["Global`*"]

T=(1/2)m x1'[t]^2+(1/2)m x2'[t]^2

V=(1/2)(3k) x1[t]^2 +(1/2)(2k) (x1[t]-x2[t])^2 //Expand

lag=T-V;

D[T,v1]

D[D[lag,x1'[t]],t] - D[lag,x1[t]]

D[D[lag,x2'[t]],t] - D[lag,x2[t]]

Tmat=m  DiagonalMatrix[{1,1}];
Tmat //MatrixForm

D[V,x1]

D[V,x2]

Vmat=k{{5,-2},{-2,2}};
Vmat //MatrixForm

mat= Vmat-Tmat ω2;
mat //MatrixForm

sol=Solve[Det[mat] ==0,ω2]

eq1=mat .{a,b} ==0 //Thread;
eq1 //Column

norm = a^2+b^2==1

eq2=Join[eq1,{norm}];
eq2 //Column

sol1=Solve[eq2 /.sol[[1]],{a,b}] //Last

ev1={a,b} /.sol1

sol2=Solve[eq2 /.sol[[2]],{a,b}]

ev2={a,b} /.Last[sol2]

eVecs={ev1,ev2} //Simplify

eVecs.eVecs

Vdiag=eVecs.Vmat.Transpose[eVecs] //Simplify ;
Vdiag//MatrixForm

Tdiag=eVecs.Tmat.Transpose[eVecs] //Simplify;
Tdiag //MatrixForm

Part B Look at motion

values={k->1,m->1};

sol

mode1=ev1 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[1]]

mode2=ev2 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[2]]

mode1 +{-1,1} /.values

Plot[mode1 +{-1,1}/.values //Re,{t,0,10}]

Plot[mode2 +{-1,1}/.values //Re,{t,0,10}]

Plot[mode1 + mode2  +{-2,2}/.values //Re,{t,0,10}]

{x1,x2}==eVecs.{n1,n2}

{x1,x2}==eVecs.{n1,n2} /.{n1->1/2,n2->1/2}

mix[a_]:=a mode1 +(1-a)mode2 /.values //Re
mix[1/2] /.{t->0} //Simplify

doPlot[a_:1/2]:=Plot[a (mode1 +{-1,1}/.values)+(1-a)(mode2 +{-1,1}/.values)
 //Re,{t,0,10}]

doPlot[1/2]

Problem #4)

Clear["Global`*"]

T=(1/2)m1 x'[t]^2+(1/2)m2 (x'[t]+L θ'[t])^2 //Expand

V=m2 g L (1-Cos[θ[t]])+ (1/2)k ((x[t])^2)//Expand

lag=T-V

rule=Cos[θ[t]]->Series[Cos[θ[t]],{θ[t],0,2}] //Normal

V=V /. rule //Simplify

D[T,x'[t]]

D[T,θ'[t]]

t11=Coefficient[ D[T,x'[t]],x'[t]]

t12=t21=Coefficient[ D[T,x'[t]],θ'[t]]

t22=Coefficient[ D[T,θ'[t]],θ'[t]]

Tmat={{t11,t12},{t21,t22}};
Tmat //MatrixForm

D[V,x[t]]

D[V,θ[t]]

v11=Coefficient[ D[V,x[t]],x[t]]

v12=v21=Coefficient[ D[V,x[t]],θ[t]]

v22=Coefficient[ D[V,θ[t]],θ[t]]

Vmat={{v11,v12},{v21,v22}};
Vmat //MatrixForm

Tmat //MatrixForm

(* x Equation *)
D[D[lag,x'[t]],t] - D[lag,x[t]] //Simplify

(* θ Equation *)
D[D[lag,θ'[t]],t] - D[lag,θ[t]]

values={m1->1, m2->1,L->1, g->1,k->2};

mat= Vmat-Tmat ω2;
mat //MatrixForm

mat=mat  //.values//Simplify;
mat //MatrixForm

sol=Solve[Det[mat] ==0,ω2] //Simplify

eq1=mat .{a,b} ==0 //Thread;
eq1 //Column

sol1=Solve[eq1 /.sol[[1]],{a,b}][[1]]   /.values //Simplify

ev1={a,b} /.sol1 /.{a->1}//Simplify

sol2=Solve[eq1 /.sol[[2]],{a,b}][[1]] /.values  //Simplify //PowerExpand

ev2={a,b} /.sol2 /.{a->1} //Simplify

eVecs=Normalize /@ {ev1,ev2} //Simplify

Vdiag=eVecs.Vmat.Transpose[eVecs]  /.values//FullSimplify ;
Vdiag//MatrixForm

Tdiag=eVecs.Tmat.Transpose[eVecs] /.values //FullSimplify ;
Tdiag //MatrixForm

Part B Look at motion

sol /.values

mode1=ev1 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[1]]  /.values

mode2=ev2 Exp[I ω t] /.{ω->Sqrt[ω2]} /.sol[[2]]/.values

mode1 +{-1,1} /.values

Plot[mode1 +{-1,1}/.values //Re,{t,0,10}]

Plot[mode2 +{-2,2}/.values //Re,{t,0,10}]

Plot[mode1 + mode2  +{-2,2}/.values //Re,{t,0,10}]

{x1,x2}==eVecs.{n1,n2}

{x1,x2}==eVecs.{n1,n2} /.{n1->1/2,n2->1/2}

mix[a_]:=a mode1 +(1-a)mode2 /.values //Re
mix[1/2] /.{t->0} //Simplify

doPlot[a_:1/2]:=Plot[a (mode1 +{-1,1}/.values)+(1-a)(mode2 +{-1,1}/.values)
 //Re,{t,0,10}]

doPlot[1/2]

Problem #5 & 6)

Clear["Global`*"]

mass= m DiagonalMatrix[{1,1}];
kmat=k {{2,-1},{-1,2}};
mass //MatrixForm
kmat //MatrixForm

mat=kmat- w2 mass;
mat //MatrixForm

Det[mat]

sol=Solve[Det[mat]==0,w2]

vec={a,b};
eqs= mat.vec ==0/.sol[[1]] //Thread

Solve[eqs,{a,b}]

ev1={1,1}

vec={a,b};
eqs= mat.vec ==0/.sol[[2]] //Thread

Solve[eqs,{a,b}]

ev2={1,-1}

Part B Look at motion

tmp1=a1  Cos[w1 t - δ1]

tmp1 //TrigExpand

tmp2= b1 Cos[w1 t] + c1 Sin[w1 t]

list1=CoefficientList[tmp2, {Cos[t w1], Sin[t w1]}]

list2=CoefficientList[tmp1 //TrigExpand, {Cos[t w1], Sin[t w1]}]

Part B Look at motion

values={k->1,m->1};

sol /.values

mode1=ev1 Exp[I w t] /.{w->Sqrt[w2]} /.sol[[1]]  /.values

mode2=ev2 Exp[I w t] /.{w->Sqrt[w2]} /.sol[[2]]  /.values

mode1 +{-1,1} /.values

Plot[mode1 +{-2,2}/.values //Re,{t,0,10}]

Plot[mode2 +{-2,2}/.values //Re,{t,0,10}]

Plot[mode1 + mode2  +{-2,2}/.values //Re,{t,0,10}]

Part B Look at motion

mode1=ev1 (b1 Cos[w t] + c1 Sin[w t])/.{w->Sqrt[w2]} /.sol[[1]]  /.values;
mode1  //MatrixForm

mode2=ev2 (b2 Cos[w t] + c2 Sin[w t])/.{w->Sqrt[w2]} /.sol[[2]]  /.values;
mode2  //MatrixForm

x[t_]=mode1 + mode2

v[t_]=D[x[t],t]

x[0]

v[0]

Find x={a,a} v={0,0}

eq1= x[0] =={a,a} //Thread

solb=Solve[eq1,{b1,b2}][[1]]

eq1v= v[0] =={0,0} //Thread

solc=Solve[eq1v,{c1,c2}][[1]]

x[t] /.solb /.solc /.{a->1}

p1=Plot[ x[t][[1]] /.solb /.solc /.{a->1},{t,0,10}];
p2=Plot[ x[t][[2]] /.solb /.solc /.{a->1},{t,0,10}];
GraphicsGrid[{{p1,p2}}]

Find x={a,0} v={0,0}

eq1= x[0] =={a,0} //Thread

solb=Solve[eq1,{b1,b2}][[1]]

eq1v= v[0] =={0,0} //Thread

solc=Solve[eq1v,{c1,c2}][[1]]

x[t] /.solb /.solc /.{a->1}

p1=Plot[ x[t][[1]] /.solb /.solc /.{a->1},{t,0,10}];
p2=Plot[ x[t][[2]] /.solb /.solc /.{a->1},{t,0,10}];
GraphicsGrid[{{p1,p2}}]

Problem 6)

Find x={0,0} v={a,a}

eq1= x[0] =={0,0} //Thread

solb=Solve[eq1,{b1,b2}][[1]]

eq1v= v[0] =={a,a} //Thread

solc=Solve[eq1v,{c1,c2}][[1]]

x[t] /.solb /.solc /.{a->1}

p1=Plot[ x[t][[1]] /.solb /.solc /.{a->1},{t,0,10}];
p2=Plot[ x[t][[2]] /.solb /.solc /.{a->1},{t,0,10}];
GraphicsGrid[{{p1,p2}}]

Find x={0,0} v={0,a}

eq1= x[0] =={0,0} //Thread

solb=Solve[eq1,{b1,b2}][[1]]

eq1v= v[0] =={0,a} //Thread

solc=Solve[eq1v,{c1,c2}][[1]]

x[t] /.solb /.solc /.{a->1}

p1=Plot[ x[t][[1]] /.solb /.solc /.{a->1},{t,0,10}];
p2=Plot[ x[t][[2]] /.solb /.solc /.{a->1},{t,0,10}];
GraphicsGrid[{{p1,p2}}]

Problem #7)

Clear["Global`*"]

x1[t_]=l1 Sin[θ1[t]];
y1[t_]=-l1 Cos[θ1[t]];

x2[t_]=l1 Sin[θ1[t]]+l2 Sin[θ2[t]];
y2[t_]=-l1 Cos[θ1[t]]-l2 Cos[θ2[t]];

x3[t_]=l1 Sin[θ1[t]]+l2 Sin[θ2[t]]+l3 Sin[θ3[t]];
y3[t_]=-l1 Cos[θ1[t]]-l2 Cos[θ2[t]]-l3 Cos[θ3[t]];

v1[t_]=Sqrt[((x1'[t])^2)+((y1'[t])^2)] //Simplify //PowerExpand

v2[t_]=Sqrt[((x2'[t])^2)+((y2'[t])^2)] //Simplify //PowerExpand

v3[t_]=Sqrt[((x3'[t])^2)+((y3'[t])^2)] //Simplify //PowerExpand

Look at single pendulum case

tKinetic1= (1/2)m1 ((v1[t])^2)//Simplify

vPot1= m1 g y1[t]

lag1=tKinetic1- vPot1;

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

Series[eq1,{θ1[t],0,1}] //Normal //Simplify

Look at double pendulum case

list={θ1[t],θ1'[t],θ2[t],θ2'[t]};
rule2=(Outer[Times,list,list]  //Flatten) ->0 //Thread

tKinetic2= (1/2)m1 ((v1[t])^2)+ (1/2)m2 ((v2[t])^2) //Simplify

vPot2= m1 g y1[t] + m2 g y2[t]

lag2=tKinetic2- vPot2;

D[lag2,θ1[t]]//Simplify

D[lag2,θ1'[t]] //Simplify

eq21=D[D[lag2,θ1'[t]],t] - D[lag2,θ1[t]]  //Simplify

eq21a=Series[eq21,{θ1[t],0,1},{θ2[t],0,1}] //Normal //ExpandAll

eq21a=eq21a   //.rule2 //Simplify

eq22=D[D[lag2,θ2'[t]],t] - D[lag2,θ2[t]]  //Simplify

eq22a=Series[eq22,{θ1[t],0,1},{θ2[t],0,1}] //Normal //ExpandAll

eq22a=eq22a   //.rule2 //Simplify

Double Pendulum: find eigensystem

(({eq21a,eq22a})/(L m)) //.{m1->m,m2->m,l1->L,l2->L} //Simplify //TableForm

mMat=L {{2,1},{1,1}};
mMat //MatrixForm

vMat=g {{2,0},{0,1}};
vMat //MatrixForm

vMat - w2 mMat //MatrixForm

sol=Solve[Det[vMat - w2 mMat ]==0,w2] //Simplify

sol1=Solve[(vMat - w2 mMat).{a,b}==0  /.sol[[1]]//Thread,{a,b}][[1]]

ev1={a,b} /.sol1 /.{a->1}   //Simplify

sol2=Solve[(vMat - w2 mMat).{a,b}==0  /.sol[[2]]//Thread,{a,b}][[1]]

ev2={a,b} /.sol2 /.{a->1}   //Simplify

evMat={ev1,ev2};
evMat.mMat.Transpose[evMat] //Simplify //MatrixForm

evMat.vMat.Transpose[evMat] //Simplify //MatrixForm

Look at triple pendulum case

list={θ1[t],θ1'[t],θ2[t],θ2'[t],θ3[t],θ3'[t]};
rule3=(Outer[Times,list,list]  //Flatten) ->0 //Thread

tKinetic3= (1/2)m1 ((v1[t])^2)+ (1/2)m2 ((v2[t])^2)+ (1/2)m3 ((v3[t])^2) //Simplify

vPot3= m1 g y1[t] + m2 g y2[t]+ m3 g y3[t]

lag3=tKinetic3- vPot3;

D[lag3,θ1[t]]//Simplify

D[lag3,θ1'[t]] //Simplify

eq31=D[D[lag3,θ1'[t]],t] - D[lag3,θ1[t]]  //Simplify

eq31a=Series[eq31,{θ1[t],0,1},{θ2[t],0,1},{θ3[t],0,1}] //Normal //ExpandAll

eq31a=eq31a   //.rule3 //Simplify

eq32=D[D[lag3,θ2'[t]],t] - D[lag3,θ2[t]]  //Simplify

eq32a=Series[eq32,{θ1[t],0,1},{θ2[t],0,1},{θ3[t],0,1}] //Normal //ExpandAll

eq32a=eq32a   //.rule3 //Simplify

eq33=D[D[lag3,θ3'[t]],t] - D[lag3,θ3[t]]  //Simplify

eq33a=Series[eq33,{θ1[t],0,1},{θ2[t],0,1},{θ3[t],0,1}] //Normal //ExpandAll

eq33a=eq33a   //.rule3 //Simplify

{eq31a,eq32a,eq33a} //TableForm

Look at double pendulum case in matrix form

eq21a
eq22a

mMat={{D[eq21a,θ1''[t]],D[eq21a,θ2''[t]]},
{D[eq22a,θ1''[t]],D[eq22a,θ2''[t]]}};
mMat //MatrixForm

kMat={{D[eq21a,θ1[t]],D[eq21a,θ2[t]]},
{D[eq22a,θ1[t]],D[eq22a,θ2[t]]}};
kMat //MatrixForm

values={m1->m,m2->m,g->1,l1->len,l2->len,len->1};

mMat //.values //MatrixForm
kMat //.values //MatrixForm

eq=Det[kMat - w2 mMat]==0 //.values

sol=Solve[eq,w2]

vec={a,b};
eq1=(kMat - w2 mMat).vec ==0 /.sol[[1]]  //.values

sol1=Solve[eq1,{a,b}] //Simplify

ev1={a,b} /.sol1[[1]] /.{a->1}

vec={a,b};
eq2=(kMat - w2 mMat).vec ==0 /.sol[[2]]  //.values

sol2=Solve[eq2,{a,b}] //Simplify

ev2={a,b} /.sol2[[1]] /.{a->1}

ev1 //N
ev2 //N