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 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 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 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 at 14:12
1  
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 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}];
ListAnimate[frames]

animation

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, 0, tMax, tMax/nFrames}];
ListAnimate[frames]

animations 2

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

share|improve this answer
    
Thank you very very much! –  Thison May 18 at 16:45

Your Answer

 
discard

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.