npendula = 4; (*In principle the code can run with any number of pendula, but some minor adjustments might be needed here and there to be sure to select the right frequencies etc.*)
\[Theta] = ToExpression[ StringJoin["\[Theta]", #] & /@ Evaluate[ToString /@ Range[npendula]]];
l = ToExpression[StringJoin["l", #] & /@ Evaluate[ToString /@ Range[npendula]]];
 m = ToExpression[StringJoin["m", #] & /@ Evaluate[ToString /@ Range[npendula]]];
 g =.;
p = {0, 0};
pos = Reap[For[j = 1, j <= npendula, j++,
    p = p + l[[j]] {Sin[\[Theta][[j]][t]], -Cos[\[Theta][[j]][t]]};
    Sow[p];
    ]][[2, 1]];
vel = D[#, t] & /@ pos;
T = FullSimplify[Sum[m[[j]]/2 (vel[[j, 1]]^2 + vel[[j, 2]]^2), {j, 1, npendula}] ];
V = g Sum[m[[j]] *pos[[j, 2]], {j, 1, npendula}];
L = T - V; (*Lagrangian of the system*)
eq = Flatten[Table[
    FullSimplify[
     (D[D[L, Evaluate[D[\[Theta][[j]][t], t]] ], t] - 
        D[L, \[Theta][[j]][t] ]) == 0
     ]
    , {j, 1, npendula}] ] /. {Join[{ g -> 1}, 
    Table[l[[j]] -> 1, {j, 1, npendula}], 
    Table[m[[j]] -> 1, {j, 1, npendula}]]}; (*Equations of motion*)
L1 = (Normal@Series[(L /. Flatten@Join[
         Table[{\[Theta][[j]][t] -> \[Epsilon] \[Theta][[j]][t]}, {j, 
           1, npendula}]
         ,
         Table[{\[Theta][[j]]'[
             t] -> \[Epsilon] \[Theta][[j]]'[t]}, {j, 1, npendula}]
         ]), {\[Epsilon], 0, 2}]) /. {\[Epsilon] -> 1}; (*Linearized Lagrangian*)
eq1 = Table[
   FullSimplify[
    (D[D[L1, \[Theta][[j]]'[t] ], t] - D[L1, \[Theta][[j]][t] ]) == 0
    ]
   , {j, 1, npendula}] // Flatten; (*Linearized equations of motion*)
f\[Theta] = 
  ToExpression[
   StringJoin["f\[Theta]", #] & /@ 
    Evaluate[ToString /@ Range[npendula]]]; (*Dummy variables*)
M = Normal@CoefficientArrays[Table[
      eq1[[j]] /. 
       Join[Table[\[Theta][[j]][t] -> f\[Theta][[j]], {j, 1, 
          npendula}],
        Table[\[Theta][[j]]'[t] -> I \[Omega] f\[Theta][[j]], {j, 1, 
          npendula}],
        Table[\[Theta][[j]]''[t] -> - \[Omega]^2 f\[Theta][[j]], {j, 
          1, npendula}] ]
      , {j, 1, npendula}], f\[Theta]][[2]];
naturalfreq1 = 
 N@Solve[Det[(M /. 
       Join[{g -> 1}, Table[l[[j]] -> 1, {j, 1 npendula}], 
        Table[m[[j]] -> 1, {j, 1 npendula}]])] == 0, \[Omega]]; (*Natural frequencies. Selecting the correct ones is probably the biggest thing to check if you use a odd number of pendula.*)
e = ToExpression[ StringJoin["e", #] & /@ Evaluate[ToString /@ Range[npendula]]]; (*Dummy variables*)
modes = Table[
  FullSimplify@ Solve[(M /. Join[{\[Omega] -> naturalfreq1[[2*j, 1, 2]], g -> 1}, Table[l[[j]] -> 1, {j, 1 npendula}], Table[m[[j]] -> 1, {j, 1 npendula}]]) . e == Table[0, {j, 1, npendula}], e]
  , {j, 1, npendula}];
orthogonalmodes = Simplify@Table[e/e1 /. modes[[j, 1]], {j, 1, npendula}]
metric = Normal@CoefficientArrays[eq1, Table[\[Theta][[j]]''[t], {j, 1, npendula}]][[2]]; (*If you are changing the number of pendula, make sure you are selecting the correct matrix.*)
metric1 = (metric /. {Join[{ g -> 1},  Table[l[[j]] -> 1, {j, 1, npendula}], Table[m[[j]] -> 1, {j, 1, npendula}]]})[[1]];
eq = Flatten[Table[
     FullSimplify[
      (D[D[L, Evaluate[D[\[Theta][[j]][t], t]] ], t] - 
         D[L, \[Theta][[j]][t] ]) == 0
      ]
     , {j, 1, npendula}] ] /. {Join[{ g -> 1}, 
     Table[l[[j]] -> 1, {j, 1, npendula}], 
     Table[m[[j]] -> 1, {j, 1, npendula}]]};
eqbound = (Join[eq, 
     Table[\[Theta][[j]][0] == 1.*orthogonalmodes[[1, j]], {j, 1, 
       npendula}], 
     Table[\[Theta][[j]]'[0] == 0, {j, 1, npendula}]]) /. {Join[{ 
      g -> 1}, Table[l[[j]] -> 1, {j, 1, npendula}], 
     Table[m[[j]] -> 1, {j, 1, npendula}]]};
vars = Table[\[Theta][[j]][t], {j, 1, npendula}];
tmax = 150;
sol = NDSolve[eqbound, vars, {t, 0, tmax}, Method -> {"EquationSimplification" -> "Residual"}] (*Solve the equations of motion.*)
solpos = (pos /. Table[l[[j]] -> 1, {j, 1, npendula}]) /. sol;
frames = Table[
   modepos = 
    Table[(pos /. Table[l[[j]] -> 1, {j, 1, npendula}]) /. 
      Table[\[Theta][[j]][t] -> orthogonalmodes[[k, j]]*
         Simplify[
           orthogonalmodes[[k]] . metric1 . 
            Evaluate[((Table[\[Theta][[j]][t], {j, 1, npendula}] /. 
                  sol) /. {t -> \[Tau]})[[1]] ] ]/(orthogonalmodes[[
             k]] . metric1 . orthogonalmodes[[k]])
       , {j, 1, npendula}], {k, 1, npendula}];
   GraphicsGrid[{{
      Graphics[{
        Line[Join[{{0, 0}}, Table[solpos[[1, j]], {j, 1, npendula}]] ],
        Disk[{0, 0}, 0.075], 
        Table[Disk[solpos[[1, j]], 0.1], {j, 1, npendula}],
        Text[
         Style["Quadruple pendulum", Bold, FontSize -> 14], {0, 3.5}]
        }, 
       PlotRange -> 
        1.1 {{-npendula, npendula}, {-npendula, npendula}}
       ]
      , SpanFromLeft,
      Graphics[{
        Line[
         Join[{{0, 0}}, Table[modepos[[1, j]], {j, 1, npendula}]] ],
        Disk[{0, 0}, 0.075], 
        Table[Disk[modepos[[1, j]], 0.1], {j, 1, npendula}],
        Text[
         Style["\!\(\*SuperscriptBox[\(1\), \(st\)]\) mode", Bold, 
          FontSize -> 14], {0, 3.5}]
        }, 
       PlotRange -> 
        1.1 {{-npendula, npendula}, {-npendula, npendula}}
       ]
      ,
      Graphics[{
        Line[
         Join[{{0, 0}}, Table[modepos[[2, j]], {j, 1, npendula}]] ],
        Disk[{0, 0}, 0.075], 
        Table[Disk[modepos[[2, j]], 0.1], {j, 1, npendula}],
        Text[
         Style["\!\(\*SuperscriptBox[\(2\), \(nd\)]\) mode", Bold, 
          FontSize -> 14], {0, 3.5}]
        }, 
       PlotRange -> 
        1.1 {{-npendula, npendula}, {-npendula, npendula}}
       ]
      }, {SpanFromAbove, SpanFromBoth,
      Graphics[{
        Line[
         Join[{{0, 0}}, Table[modepos[[3, j]], {j, 1, npendula}]] ],
        Disk[{0, 0}, 0.075], 
        Table[Disk[modepos[[3, j]], 0.1], {j, 1, npendula}],
        Text[
         Style["\!\(\*SuperscriptBox[\(3\), \(rd\)]\) mode", Bold, 
          FontSize -> 14], {0, 3.5}]
        }, 
       PlotRange -> 
        1.1 {{-npendula, npendula}, {-npendula, npendula}}
       ]
      ,
      Graphics[{
        Line[
         Join[{{0, 0}}, Table[modepos[[4, j]], {j, 1, npendula}]] ],
        Disk[{0, 0}, 0.075], 
        Table[Disk[modepos[[4, j]], 0.1], {j, 1, npendula}],
        Text[
         Style["\!\(\*SuperscriptBox[\(4\), \(th\)]\) mode", Bold, 
          FontSize -> 14], {0, 3.5}]
        }, 
       PlotRange -> 
        1.1 {{-npendula, npendula}, {-npendula, npendula}}
       ]
      }}, Frame -> All, ImageSize -> 600]
   , {t, 0, tmax/1, 0.3}];
ListAnimate[frames] (*Plot everything.*)