File:Action-Angle variables for a coupled harmonic oscillator.gif

m = 1; k = 1/2; k1 = 2; \[Omega]1 = Sqrt[(2 k)/m]; \[Omega]2 = Sqrt[(2 k1)/m];
H1[t_] := p1[t]^2/(2 m) + k y1[t]^2; H2[t_] := p2[t]^2/(2 m) + k1 y2[t]^2; (*A coupled harmonic oscillator can be rewritten as two uncoupled harmonic oscillators via a change of variables*)
sol1 = FullSimplify[ DSolve[{p1'[t] == D[H1[t], y1[t] ], y1'[t] == -D[H1[t], p1[t]]}, {p1[t], y1[t]}, t] ][[1]];
sol2 = FullSimplify[ DSolve[{p2'[t] == D[H2[t], y2[t] ], y2'[t] == -D[H2[t], p2[t]]}, {p2[t], y2[t]}, t] ][[1]];
ctot[{\[Theta]_, \[Phi]_}] := {(2 + Cos[\[Theta]]) Cos[\[Phi]], (2 + Cos[\[Theta]]) Sin[\[Phi]], Sin[\[Theta]]}; (*Equations of a torus*)
spring[p1_, p2_, ncoils_, leads_, amplitude_] := (
   p22 = p2 - p1;
   \[Theta] = If[Norm[p22] > 0, ArcTan[p22[[2]], p22[[1]]] - \[Pi]/2, 0];
   p23 = Norm[p22] {1, 0};
   ((# + p1) & /@ RotationTransform[-\[Theta]][ If[ Norm[p2 - p1] > 2 leads, Join[{{0, 0}}, {{0, 0} + {leads, 0}}, Table[{(Norm[p23] - 2 leads)/ncoils j, amplitude (-1)^j} + {0,0} + {leads, 0}, {j, 1, ncoils - 1}], {p23 - {leads, 0}}, {p23}] 
       ,
       Join[{{0, 0}, p23/2, p23/2 + {0, amplitude}, p23/2 - {0, amplitude}, p23/2, p23}]
       ]
      ]) );
x1[t_] := (y1[t] + y2[t])/Sqrt[2] /. sol1 /. sol2 /. {C[1] -> 0, C[2] -> 1}; x2[t_] := (y1[t] - y2[t])/Sqrt[2] /. sol1 /. sol2 /. {C[1] -> 0, C[2] -> 1}; (*Go back to the original set of variables*)
frames = Table[
   Grid[{{
      Graphics[{
        Black, PointSize[0.02], Point[{{0, 0}, {3 L, 0}}],
        Line[spring[{0, 0}, {x1[t] + L, 0}, 20, 0.2, 0.1]],
        Line[spring[{x2[t] + 2 L, 0}, {x1[t] + L, 0}, 20, 0.2, 0.1]],
        Line[spring[{x2[t] + 2 L, 0}, {3 L, 0}, 20, 0.2, 0.1]],
        Style[
         Text["\!\(\*SubscriptBox[\(x\), \(1\)]\)", {x1[t] + 
            L, -0.5}], Black, Bold, Large],
        Style[
         Text["\!\(\*SubscriptBox[\(x\), \(2\)]\)", {x2[t] + 
            2 L, -0.5}], Black, Bold, Large],
        Red, PointSize[0.03], Point[{{x1[t] + L, 0}, {x2[t] + 2 L, 0}}]
        }, ImageSize -> Large]
      }, {
      Show[
       ParametricPlot3D[
        ctot[{\[Theta], \[Phi]}], {\[Theta], 0, 2 \[Pi]}, {\[Phi], 0, 
         2 \[Pi]}, Mesh -> None, Boxed -> False, Ticks -> None, 
        AxesLabel -> \
{"(2+\!\(\*SqrtBox[FractionBox[SuperscriptBox[\((\*SubscriptBox[\
OverscriptBox[\(x\), \(.\)], \(1\)] + \
\*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(2\)])\), \(2\)], \
\(\*SuperscriptBox[\((\*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(1\
\)] + \*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(2\)])\), \(2\)] + \
\*SuperscriptBox[\((\*SubscriptBox[\(x\), \(1\)] + \*SubscriptBox[\(x\
\), \(2\)])\), \(2\)]\\\ \*SuperscriptBox[SubscriptBox[\(\[Omega]\), \
\(1\)], \(2\)]\)]]\))\!\(\*SqrtBox[FractionBox[SuperscriptBox[\((\*\
SubscriptBox[OverscriptBox[\(x\), \(.\)], \(1\)] - \
\*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(2\)])\), \(2\)], \
\(\*SuperscriptBox[\((\*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(1\
\)] - \*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(2\)])\), \(2\)] + \
\*SuperscriptBox[\((\*SubscriptBox[\(x\), \(1\)] - \*SubscriptBox[\(x\
\), \(2\)])\), \(2\)]\\\ \*SuperscriptBox[SubscriptBox[\(\[Omega]\), \
\(2\)], \(2\)]\)]]\)", 
          "(2+\!\(\*SqrtBox[FractionBox[SuperscriptBox[\((\*\
SubscriptBox[OverscriptBox[\(x\), \(.\)], \(1\)] + \
\*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(2\)])\), \(2\)], \
\(\*SuperscriptBox[\((\*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(1\
\)] + \*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(2\)])\), \(2\)] + \
\*SuperscriptBox[\((\*SubscriptBox[\(x\), \(1\)] + \*SubscriptBox[\(x\
\), \(2\)])\), \(2\)]\\\ \*SuperscriptBox[SubscriptBox[\(\[Omega]\), \
\(1\)], \(2\)]\)]]\))\!\(\*SqrtBox[FractionBox[\(\*SuperscriptBox[\
SubscriptBox[\(\[Omega]\), \(2\)], \(2\)] \
\*SuperscriptBox[\((\*SubscriptBox[\(x\), \(1\)] - \*SubscriptBox[\(x\
\), \(2\)])\), \(2\)]\), \
\(\*SuperscriptBox[\((\*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(1\
\)] - \*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(2\)])\), \(2\)] + \
\*SuperscriptBox[\((\*SubscriptBox[\(x\), \(1\)] - \*SubscriptBox[\(x\
\), \(2\)])\), \(2\)]\\\ \*SuperscriptBox[SubscriptBox[\(\[Omega]\), \
\(2\)], \(2\)]\)]]\)", 
          "\!\(\*SqrtBox[FractionBox[\(\*SuperscriptBox[SubscriptBox[\
\(\[Omega]\), \(1\)], \(2\)] \
\*SuperscriptBox[\((\*SubscriptBox[\(x\), \(1\)] + \*SubscriptBox[\(x\
\), \(2\)])\), \(2\)]\), \
\(\*SuperscriptBox[\((\*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(1\
\)] + \*SubscriptBox[OverscriptBox[\(x\), \(.\)], \(2\)])\), \(2\)] + \
\*SuperscriptBox[\((\*SubscriptBox[\(x\), \(1\)] + \*SubscriptBox[\(x\
\), \(2\)])\), \(2\)]\\\ \*SuperscriptBox[SubscriptBox[\(\[Omega]\), \
\(1\)], \(2\)]\)]]\)"}, PlotStyle -> {Gray, Opacity[0.25]}, 
        Lighting -> "Neutral", LabelStyle -> {Tiny, Black, Bold}, 
        ImagePadding -> All, ImageSize -> Large, AxesStyle -> Black]
       ,
       ParametricPlot3D[
        ctot[{\[Omega]1 \[Tau], \[Omega]2 \[Tau]}] , {\[Tau], t - 5, 
         t}, ColorFunction -> 
         Function[{x, y, z, \[Tau]}, 
          Directive[ColorData["SunsetColors"][\[Tau]/\[Pi]] , 
           Opacity[\[Tau]^3] ]] ]
       ]
      }}]
   , {t, 0, 2 \[Pi] - \[Pi]/10, \[Pi]/10}];
ListAnimate[frames]