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.

I wrote a little program to use Newton's Law of Universal Gravitation to animate 3 planets orbiting a central star, but I have run into a problem. Here is the code that I used to create the program (I apologise about the messyness):

orbit1 = NDSolve[{
x''[t] == (-(6.672*10^-11) (7*10^17) x[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
y''[t] == (-(6.672*10^-11) (7*10^17) y[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
z''[t] == (-(6.672*10^-11) (7*10^17) z[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
x[0] == 1000, y[0] == 1000, z[0] == 1000, 
x'[0] == 0, y'[0] == -100, z'[0] == 0}, {x[t], y[t], z[t]}, {t, 0, 1000}];

orbit2 = NDSolve[{
x''[t] == (-(6.672*10^-11) (7*10^17) x[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
y''[t] == (-(6.672*10^-11) (7*10^17) y[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
z''[t] == (-(6.672*10^-11) (7*10^17) z[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
x[0] == 500, y[0] == -1000, z[0] == -1000, 
x'[0] == -110, y'[0] == 100, z'[0] == 0}, {x[t], y[t], z[t]}, {t, 0, 1000}];

orbit3 = NDSolve[{
x''[t] == (-(6.672*10^-11) (7*10^17) x[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
y''[t] == (-(6.672*10^-11) (7*10^17) y[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
z''[t] == (-(6.672*10^-11) (7*10^17) z[t])/(x[t]^2 + y[t]^2 + z[t]^2)^(3/2), 
x[0] == 0, y[0] == 100, z[0] == 500, 
x'[0] == 350, y'[0] == -100, z'[0] == 0}, {x[t], y[t], z[t]}, {t, 0, 1000}];

orbitplotunion = 
 Animate[Show[{ParametricPlot3D[{{x[t], y[t], z[t]} /. 
   orbit1, {x[t], y[t], z[t]} /. orbit2, {x[t], y[t], z[t]} /. 
   orbit3}, {t, 0, a}, 
 PlotRange -> {{-1600, 1600}, {-1600, 1600}, {-1600, 1600}}, 
 AxesLabel -> {x, y, z}], 
Graphics3D[{Yellow, Sphere[{0, 0, 0}, 100], Green, 
  Sphere[{x[t], y[t], z[t]} /. orbit1 /. t -> a, 50], Blue, 
  Sphere[{x[t], y[t], z[t]} /. orbit2 /. t -> a, 50], Purple, 
  Sphere[{x[t], y[t], z[t]} /. orbit3 /. t -> a, 50]}]}], {a, 0, 
Infinity, 0.1}]

As you can see, there is one orbit calculation for each planet, and these are then animated.

Now, my first problem has to do with NDSolve[] and its calculation from t = 0 to t = 1000, which means that the animation will break once t = 1000 is hit.

Is there a way to allow the animation to go on indefinitely, instead of having to reset it once t = 1000 comes along?

Secondly, the planets and their trailing lines start out orbiting the star perfectly in the beginning, but over time, the trailing lines of each planet become more and more jagged, until eventually the animation looks horrible.

Here is what the orbits look like in the beginning:

orbits at start

And this is what they look like some time later:

later orbits

If anyone knows how to solve the jagged line problem (maybe by editing the amount of trailing line that is allowed, so that it doesn't redraw over itself continually) and knows how to make the animation continue indefinitely, I would love to know.

Regards,

Alex


EDIT:

After looking at the $n$-body wiki page, I thought I'd give it a go and start out small with a simple Earth-Sun simulation, and if that worked, then move my way up to $3,4,5,\dots$ bodies as well. Unfortunately, as expected, I seem to have run into a problem almost immediately. Here is the code that I am currently using:

G = 6.672*10^-11;
m1 = 5.972*10^24; (* mass of Earth *)
m2 = 1.989*10^30; (* mass of Sun *)

orbitearthsun = NDSolve[{
x1''[t] == -((G (m2) (x1[t] - x2[t]))/Abs[x1[t] - x2[t]]^3),
y1''[t] == -((G (m2) (y1[t] - y2[t]))/Abs[y1[t] - y2[t]]^3),
x2''[t] == -((G (m1) (x2[t] - x1[t]))/Abs[x2[t] - x1[t]]^3),
y2''[t] == -((G (m1) (y2[t] - y1[t]))/Abs[y2[t] - y1[t]]^3),
x1[0] == 1000, x2[0] == 0, x1'[0] == 100, x2'[0] == 0,
y1[0] == 1000, y2[0] == 0, y1'[0] == 100, y2'[0] == 0},
{x1[t], x2[t], y1[t], y2[t]}, {t, 0, 1000}]

NDSolve::ndsz: At t == 3.049010336028579`*^-6, step size is effectively zero;
singularity or stiff system suspected. >>

Does this occur because the denominator of each term becomes zero as they collide?

share|improve this question
1  
How about adding PerformanceGoal -> "Quality"... –  cormullion May 11 '13 at 12:40
1  
Very cool looking! I wonder how you might edit it to get your planets to interact with one another? –  Mark McClure May 11 '13 at 13:22
1  
Seems like it'll be an extremely tricky undertaking Mark. en.wikipedia.org/wiki/N-body_problem –  user7388 May 11 '13 at 13:31
3  
The jagged lines are because ParametricPlot3D isn't sampling the function finely enough when you go to large values of a. Try using something like {t, Max[0, a-50], a} as the parameter range for the plot. The 50 determines the length of the trail. –  Simon Woods May 11 '13 at 13:50
1  
Roger that, I've since made some nice progress on the n-body problem but will ask if I run into anymore troubles. Thanks for your help guys, I appreciate it. –  user7388 May 13 '13 at 8:12

2 Answers 2

up vote 15 down vote accepted

Some frames from my version of the animation:

orbiting planets

Here's the code I used:

orbit[posStart_?VectorQ, derStart_?VectorQ] := 
     Block[{c = -Rationalize[6.672*^-11*7*^17], x, y, z, t},
           {x, y, z} /. First @ NDSolve[
           Join[Thread[{x''[t], y''[t], z''[t]} == 
                       c {x[t], y[t], z[t]}/Norm[{x[t], y[t], z[t]}]^3], 
                Thread[{x[0], y[0], z[0]} == posStart], 
                Thread[{x'[0], y'[0], z'[0]} == derStart]],
           {x, y, z}, {t, 0, ∞},
           Method -> {"EventLocator", Direction -> 1,
                      "Event" -> {x'[t], y'[t], z'[t]}.({x[t], y[t], z[t]} - posStart),
                      "EventAction" :> Throw[Null, "StopIntegration"], 
                      Method -> {"SymplecticPartitionedRungeKutta", 
                                 "PositionVariables" -> {x[t], y[t], z[t]}}}, 
           WorkingPrecision -> 20]]

{x[1], y[1], z[1]} = orbit[{1000, 1000, 1000}, {0, -100, 0}];
tf1 = x[1]["Domain"][[1, -1]];

{x[2], y[2], z[2]} = orbit[{500, -1000, -1000}, {-110, 100, 0}];
tf2 = x[2]["Domain"][[1, -1]];

{x[3], y[3], z[3]} = orbit[{0, 100, 500}, {350, -100, 0}];
tf3 = x[3]["Domain"][[1, -1]];

orbit1[t_] := Through[{x[1], y[1], z[1]}[tf1 Mod[t/tf1, 1]]];
orbit2[t_] := Through[{x[2], y[2], z[2]}[tf2 Mod[t/tf2, 1]]];
orbit3[t_] := Through[{x[3], y[3], z[3]}[tf3 Mod[t/tf3, 1]]];

Animate[Show[
        ParametricPlot3D[{orbit1[t], orbit2[t], orbit3[t]}, {t, -$MachineEpsilon, a}], 
        Graphics3D[{{Yellow, Sphere[{0, 0, 0}, 100]},
                    {Green, Sphere[orbit1[a], 50]},
                    {Blue, Sphere[orbit2[a], 50]},
                    {Purple, Sphere[orbit3[a], 50]}}],
                   AxesLabel -> {x, y, z}, 
                   PlotRange -> {{-1600, 1600}, {-1600, 1600}, {-1600, 1600}}],
        {a, 0, ∞, 2}]

The only non-basic idea here is the use of event detection in conjunction with a symplectic integrator; briefly, one should certainly use a symplectic integrator when integrating a Hamiltonian system, and that you can use event detection to detect periodic behavior when the solutions of a system of differential equations is known to be periodic.

share|improve this answer
    
(For large values of t, the curves might start looking jittery due to roundoff; cranking up WorkingPrecision might be necessary. Simon's comment on sampling is also relevant.) –  J. M. May 11 '13 at 15:56
    
Wouldn't a symplectic integrator approach only be necessary if you're looking to model complex systems? –  user7388 May 11 '13 at 18:53
1  
@user7388: A symplectic integrator preserves the symmetry of any Hamiltonian system. –  chris May 11 '13 at 19:53
    
(If anyone wants to write a WhenEvent[] version, feel free; I don't have version 9 for doing that.) –  J. M. May 16 '13 at 11:06

This is too long for a comment, but it isn't an answer.

Perhaps you would like to consider a more compact way to write down your equations:

m = -(6.672*10^-11) (7*10^17) ;
st = {{1000, 1000, 1000, 0, -100, 0},
      {500, -1000, -1000, -110, 100, 0},
      {0, 100, 500, 350, -100, 0}};
r = {x @ t, y @ t, z @ t};

o[n_] := NDSolve[Join[{Equal @@ Join[(D[r, {t, 2}]/r), {m/Norm @ r^3}]},
                       Thread[{x[0], y[0], z[0], x'[0], y'[0], z'[0]} == st[[n]]]],
                 r, {t, 0, 1000}]

And the plotting part gets to:

d = {d1, d2, d3} = Evaluate[r /. o /@ {1, 2, 3}];
Animate[Show[{
   ParametricPlot3D[d /. t -> u, {u, 0, a}, PlotRange->1600 {{-1, 1}, {-1, 1}, {-1, 1}}], 
   Graphics3D[ MapThread[{#1, Sphere[#2 /. t -> a, #3]} &, {{Yellow, Green, Blue, Purple}, 
                                            {{0, 0, 0}, d1, d2, d3}, 50 {2, 1, 1, 1}}]]}], 
{a, 0, 100, 0.1}]
share|improve this answer
    
Thanks belisarius, I'll have a look at compacting the code. –  user7388 May 11 '13 at 14:31

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.