Small Oscillations of the n-pendulum

Author: Ryan Rubenzahl  (12/10/2016)

n=30;

g=9.81;a=1/n;

Do[θ[i][t_]=q[i][t],{i,1,n}];
Do[θd[i][t_]=q[i]'[t],{i,1,n}]
Do[θdd[i][t_]=q[i]''[t],{i,1,n}]
Do[BCs[i]=Sequence@@{θ[i][0]==π/16,θd[i][0]==0},{i,1,n}]
BC = Table[BCs[i],{i,1,n}];

Do[EquationsOfMotion[k]=Sum[(n-Max[k,j]+1)*θdd[j][t],{j,1,n}]+(n-k+1)*(g/a)θ[k][t]==0,{k,1,n}]

Table[EquationsOfMotion[i],{i,1,n}];

sol=NDSolve[{Table[EquationsOfMotion[i],{i,1,n}],BC},Table[θ[i][t],{i,1,n}],{t,0,100},Method->{"EquationSimplification"->"Residual"}]

Do[thetas[i][t_]:=Evaluate[θ[i][t]/.sol][[1]],{i,n}]
Do[xpos[j][t_]:=Evaluate[Sum[a*Sin[thetas[i][t]],{i,1,j}]],{j,1,n}]
Do[ypos[j][t_]:=Evaluate[Sum[a*Cos[thetas[i][t]],{i,1,j}]],{j,1,n}]
Do[lines[i][t_]:=Evaluate[Line[{{xpos[i][t],-ypos[i][t]},{xpos[i+1][t],-ypos[i+1][t]}}]],{i,1,n-1}]
Do[bobs[i][t_]:=Evaluate[Disk[{xpos[i][t],-ypos[i][t]},a/5]],{i,1,n}]

TotalLength=n*a;
TotalLength=TotalLength+TotalLength/10;
Animate[Plot[{0},{x,-TotalLength,TotalLength},PlotRange->{-TotalLength,0},Ticks->None,Epilog->{Line[{{0,0},{xpos[1][t],-ypos[1][t]}}],Table[lines[i][t],{i,1,n-1}],Table[bobs[i][t],{i,1,n}]}, AspectRatio->1/2],{t,0,100,AnimationRate->1,Appearance->"Open",AnimationRunning->False}]