Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}]
share|improve this question
I see one equation that has two == in it. Probably two equations that you have inadvertently merged. Ans it may be that in {phi1[t], phi2[t], t} you may have to drop the t part. –  Sjoerd C. de Vries May 18 '14 at 13:49
Try this instead: 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]] == 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}]. Since this seems to be a case of a simple error I'm going to vote to close this question (no worries, nothing negative). –  Sjoerd C. de Vries May 18 '14 at 13:53
For a great introduction to Mathematica, have a look at Leonid Shifrin's free book. An extensive overview of additional Mathematica material can be found here. –  Sjoerd C. de Vries May 18 '14 at 13:57
x1[t] := Evaluate[Sin[phi1[t]] /. sol] should be x1[t_] = Evaluate[Sin[phi1[t]] /. sol]. And y1,x2, and y2 should be similarly changed. ParametricPlot[Evaluate[{{x1[t], y1[t]}, {x2[t], y2[t]}} /. sol], {t, 0, 10}] should be ParametricPlot[{{x1[t], y1[t]}, {x2[t], y2[t]}}], {t, 0, 10}] –  m_goldberg May 18 '14 at 14:12
You'll probably get more meaningful graphical results if you use a dynamic rather than a static plot. E.g.: Manipulate[ ParametricPlot[Evaluate[{{x1[t], y1[t]}, {x2[t], y2[t]}}], {t, 0, s}, PlotRange -> {{-2, 2}, {-2, 1.5}}], {{s, 0.00001}, 0, 10}] –  murray May 18 '14 at 15:43

1 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}];


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:


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}];

animations 2

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


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.

share|improve this answer
Thank you very very much! –  Thison May 18 '14 at 16:45
@Jens : I wonder how it is possible to add the trajectory of the pendulum ( I mean the trajectory of the red disks in the animations ) into the animation.Could you please help me ? –  Shellp Nov 5 '14 at 20:31

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.