Animation of double pendulum

Question

Sadly, I am a completely newbie. I am studying Physics and in our theoretical physics class we got the task to solve the double pendulum using Mathematica. We just got the program, but no introduction how to use it. So I am trying to animate the parametric plot. Any ideas?

sol = NDSolve[{
  2*phi1''[t] + phi2''[t]*Cos[phi1[t] - phi2[t]] + 
    phi2'[t]^(2)*Sin[phi1[t] - phi2[t]] + 2*9.81*Sin[phi1[t]] == 
  phi2''[t] + phi1''[t]*Cos[phi1[t] - phi2[t]] - 
    phi1'[t]^(2)*Sin[phi1[t] - phi2[t]] + 9.81*Sin[phi2[t]] == 0, 
  phi1[0] == Pi/2, phi2[0] == Pi, phi1'[0] == 0, phi2'[0] == 0}, 
  {phi1[t], phi2[t], t}, {t, 0, 10}] 

Plot[{phi1[t], phi2[t]} /. sol, {t, 0, 10}]

x1[t] := Evaluate[Sin[phi1[t]] /. sol]
y1[t] := Evaluate[-Cos[phi1[t]] /. sol]
x2[t] := Evaluate[Sin[phi1[t]] + Sin[phi2[t]] /. sol]
y2[t] := Evaluate[-(Cos[phi1[t]] + Cos[phi2[t]]) /. sol]

ParametricPlot[Evaluate[{{x1[t], y1[t]}, {x2[t], y2[t]}} /. sol], {t, 0, 10}]

Answer 1

After correcting the syntax errors in the original code, the actual question can be addressed:

How to display the four variables x1[t]...y2[t] as an animation in a way that conveys their meaning? The basic idea is to use ListAnimate on a list of frames that I define below:

Clear[phi1, phi2, t];
sol = First[
   NDSolve[{2*phi1''[t] + phi2''[t]*Cos[phi1[t] - phi2[t]] + 
       phi2'[t]^(2)*Sin[phi1[t] - phi2[t]] + 2*9.81*Sin[phi1[t]] == 0,
      phi2''[t] + phi1''[t]*Cos[phi1[t] - phi2[t]] - 
       phi1'[t]^(2)*Sin[phi1[t] - phi2[t]] + 9.81*Sin[phi2[t]] == 0, 
     phi1[0] == Pi/2, phi2[0] == Pi, phi1'[0] == 0, 
     phi2'[0] == 0}, {phi1[t], phi2[t]}, {t, 0, 10}]];

x1[t_] := Evaluate[Sin[phi1[t]] /. sol]
y1[t_] := Evaluate[-Cos[phi1[t]] /. sol]
x2[t_] := Evaluate[Sin[phi1[t]] + Sin[phi2[t]] /. sol]
y2[t_] := Evaluate[-(Cos[phi1[t]] + Cos[phi2[t]]) /. sol]

frames = Table[
   Graphics[{Gray, Thick, 
     Line[{{0, 0}, {x1[t], y1[t]}, {x2[t], y2[t]}}], Darker[Blue], 
     Disk[{0, 0}, .1], Darker[Red], Disk[{x1[t], y1[t]}, .1], 
     Disk[{x2[t], y2[t]}, .1]}, 
    PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}],
   {t, 0, 10, .1}];
ListAnimate[frames]

The list frames is generated in a Table where the time t is incremented in steps of .1 through the interval of the NDSolve solution. The drawing depicts the two xy coordinates as red disks, connected to each other and to the origin (the pivot) by gray lines. These are the rigid rods connecting the masses.

In terms of code, I changed the definition of sol to remove one level of List so that the resulting definition of the coordinate functions doesn't end up returning lists, but numbers. These coordinate definitions should also have a t_ instead of x1[t] etc. to be consistent with the meaning of the := (SetDelayed) in their definition (alternatively one could also have done those definitions in a simpler way without :=).

Edit: more general solution

To get more out of this animation, it may be of interest to change the simulation so that it can handle different masses, and different pendulum lengths. But this means you need to understand how the differential equation must be modified. For many constrained mechanics problems, including the double pendulum, the Lagrange formalism is the most efficient way to set up the equations of motion. Mathematica has a VariationalMethods package that helps to automate most of the steps. This is extremely useful for solving mechanics problems because it lets you focus on the physics while taking over the mathematical work of getting from the physics to the equations of motion.

Here is how it works. The following few lines contain all the steps from defining the variables to establishing and then solving the equations of motion, with arbitrary parameters for the pendulum geometry. Here s1, s2 are the lengths of the two rods connecting the masses m1, m2. Most of the following equations can be done entirely in terms of the vectorial positions r1 and r2 without thinking about their components too much. The actual variables we solve for are ϕ1, ϕ2. To show the effect of doubling one of the masses and halving one of the pendulum rods, I then set those variables before calling NDSolve:

Needs["VariationalMethods`"]

Clear[s1, s2, ϕ1, ϕ2, t, g, m1, m2];
variables = {ϕ1[t], ϕ2[t]};
r1 = s1 {Sin[ϕ1[t]], -Cos[ϕ1[t]]};
r2 = r1 + s2 {Sin[ϕ2[t]], -Cos[ϕ2[t]]};

lagrangian = 
  m1/2 D[r1, t].D[r1, t] + m2/2 D[r2, t].D[r2, t] - 
   g {0, 1}.(m1 r1 + m2 r2);
eqs = EulerEquations[lagrangian, variables, t];

s1 = 1;
s2 = 1/2;
m1 = 2;
m2 = 1;
g = 9.81;
tMax = 10;
initial = {ϕ1[0] == Pi/2, ϕ2[0] == Pi, 
           ϕ1'[0] == Pi/2, ϕ2'[0] == 0};
sol = First[NDSolve[Join[eqs, initial], variables, {t, 0, tMax}]];    
nFrames = 100;
frames = Table[Graphics[{Gray, Thick,
      Line[{{0, 0}, r1, r2}],
      Darker[Blue], Disk[{0, 0}, .1],
      Darker[Red], Disk[r1, m1/10], Disk[r2, m2/10]} /. sol,
    PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}},
    Background -> Lighter[Cyan]], 
  {t, tMax/nFrames, tMax, tMax/nFrames}];
ListAnimate[frames]

Note that the differential equations we're solving are never written down by hand anywhere. They are automatically generated by EulerEquations.

Edit

As requested in the comment, the solution above can also be re-used to plot a trace of the trajectory together with the pendulum graphic:

traceFrames = 
  Map[ParametricPlot[r2 /. sol, {t, 0, #}, 
     PlotStyle -> Directive[Thick, Lighter[Purple]], 
     PlotPoints -> Floor[(2 nFrames #)/tMax], 
     PlotRange -> {{-2.5, 2.5}, {-2.5, 2.5}}, 
     Background -> Lighter[Cyan]] &, 
   Range[tMax/nFrames, tMax, tMax/nFrames]];

ListAnimate[MapThread[Show, {traceFrames, frames}]]

Here I re-use the list of frames from above and Show it together with a new list traceFrames in which the ParametricPlot of the end point r2 is shown up to the intermediate time corresponding to the instants used in frame. Since ParametricPlot uses interpolation, we have to make sure that the plots for long times don't lose details, by increasing the number of PlotPoints with time.

Answer 2

This is not a direct answer to your problem, but rather a generalization of @Jens' code from double to n-tuple pendulums. Meaning you can also use it for double pendulums if you like. I'm providing it due to popular demand.

Needs["VariationalMethods`"]

n = 10; (* number of pendulum segments *)
rate = 10; (* animation frame rate *)
Clear[s, ϕ, t, g, m];
variables = ϕ[#][t] & /@ Range[n];
r[1] = s[1] {Sin[ϕ[1][t]], -Cos[ϕ[1][t]]};
For[i = 2, i <= n, i++,
 r[i] = r[i - 1] + s[i] {Sin[ϕ[i][t]], -Cos[ϕ[i][t]]};
 ]
lagrangian = 
  Sum[m[i]/2 D[r[i], t].D[r[i], t], {i, n}] - 
   g {0, 1}.(Sum[m[i] r[i], {i, n}]);
eqs = EulerEquations[lagrangian, variables, t];

(s[#] = 1) & /@ Range[n]; (* this is a placeholder, you can still specify individual *)
(m[#] = 1) & /@ Range[n]; (* lengths and masses *)
g = 9.81;
tMax = 10;
initial = {Sequence @@ ((ϕ[#][0] == Pi/2) & /@ Range[n]), 
   Sequence @@ ((ϕ[#]'[0] == 0) & /@ Range[n])};

sol = First[
   NDSolve[Join[eqs, initial], variables, {t, 0, tMax}, 
    Method -> {"EquationSimplification" -> "Residual"}]];

nFrames = tMax*rate;
frames = Table[
   Graphics[{Gray, Thick, 
      Line[{{0, 0}, Sequence @@ (r[#] & /@ Range[n])}], Darker[Blue], 
      Disk[{0, 0}, .1], Darker[Red], 
      Sequence @@ (Disk[r[#], m[#]/10] & /@ Range[n])} /. sol, 
    PlotRange -> {{-(n + 1), n + 1}, {-(n + 1), (n + 1)/2}}, 
    Background -> Lighter[Cyan]], {t, tMax/nFrames, tMax, 
    tMax/nFrames}];
traceFrames = 
  Map[ParametricPlot[r[n] /. sol, {t, 0, #}, 
     PlotStyle -> Directive[Thick, Lighter[Purple]], 
     PlotPoints -> Floor[(2 nFrames #)/tMax], 
     PlotRange -> {{-(n + 1), n + 1}, {-(n + 1), (n + 1)/2}}, 
     Background -> Lighter[Cyan]] &, 
   Range[tMax/nFrames, tMax, tMax/nFrames]];

ListAnimate[MapThread[Show, {traceFrames, frames}], 
 AnimationRate -> rate, AnimationRunning -> False]

This is a trajectory of a pendulum with 10 segments of equal lengths and masses right at the time it completes its first swing:

Now, if you want, you can further expand the code to include perhaps drag or spring forces between the segments.

Answer 3

I just modified the code in the graphic part so that it adapts to the different values of S1 and S2, since when establishing values higher than 2.5, the pendulums were no longer visible.

Needs["VariationalMethods`"]

Clear[s1, s2, \[Phi]1, \[Phi]2, t, g, m1, m2];
variables = {\[Phi]1[t], \[Phi]2[t]};
r1 = s1 {Sin[\[Phi]1[t]], -Cos[\[Phi]1[t]]};
r2 = r1 + s2 {Sin[\[Phi]2[t]], -Cos[\[Phi]2[t]]};

lagrangian = 
  m1/2 D[r1, t].D[r1, t] + m2/2 D[r2, t].D[r2, t] - 
   g {0, 1}.(m1 r1 + m2 r2);
eqs = EulerEquations[lagrangian, variables, t];

s1 = 4;
s2 = 4;
m1 = 3;
m2 = 6;
g = 9.81;
tMax = 10;
initial = {\[Phi]1[0] == Pi/2, \[Phi]2[0] == Pi, \[Phi]1'[0] == 
    Pi/2, \[Phi]2'[0] == 0};
sol = First[NDSolve[Join[eqs, initial], variables, {t, 0, tMax}]];
nFrames = 100;
frames = Table[
   Graphics[{Gray, Thick, Line[{{0, 0}, r1, r2}], Darker[Blue], 
      Disk[{0, 0}, .1], Darker[Red], Disk[r1, m1/10], Darker[Blue], 
      Disk[r2, m2/10]} /. sol, 
    PlotRange -> {{-(s1 + s2 + 1), (s1 + s2 + 
         1)}, {-(s1 + s2 + 1), (s1 + s2 + 1)}}, 
    Background -> Lighter[Cyan]], {t, tMax/nFrames, tMax, 
    tMax/nFrames}];
traceFrames = 
  Map[ParametricPlot[r2 /. sol, {t, 0, #}, 
     PlotStyle -> Directive[Thick, Lighter[Blue]], 
     PlotPoints -> Floor[(2 nFrames #)/tMax], 
     PlotRange -> {{-(s1 + s2 + 1), (s1 + s2 + 
          1)}, {-(s1 + s2 + 1), (s1 + s2 + 1)}}, 
     Background -> Lighter[Yellow]] &, 
   Range[tMax/nFrames, tMax, tMax/nFrames]];
ListAnimate[MapThread[Show, {traceFrames, frames}]]