3D orbits and inaccuracy over time

Question

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?

Hmm, PerformanceGoal -> "Quality" does work very well, but unfortunately the program crawls to a standstill within 2 minutes because the lines are still being redraw on top of each other. Is there a way to make the trailing lines disappear after a designated amount of time or distance? — InquisitiveInquirer, May 11 '13 at 12:51
Odd. Mine's been running at the same speed for 15 minutes... The lines are slightly thicker towards the end of the run - but I think that looks more realistic :) — cormullion, May 11 '13 at 13:01
My laptop is pretty old, so that's most likely the problem :/ — InquisitiveInquirer, May 11 '13 at 13:07
It seems one needs to give mathematica a chance every now and then to re-render what as been plotted already? — chris, May 11 '13 at 13:17
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
Seems like it'll be an extremely tricky undertaking Mark. http://en.wikipedia.org/wiki/N-body_problem — InquisitiveInquirer, May 11 '13 at 13:31
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
Thanks Simon, that made the program run much faster! Now all I need to figure out is how to make the orbits work to t=Infinity without messing up. — InquisitiveInquirer, May 11 '13 at 14:33
J.M.'s answer is very nice I think if you are ignoring the interactions between the masses orbiting the central mass. What would be interesting to see is how they do influence each other over time... — SEngstrom, May 11 '13 at 16:22
The two-body problem is a nice one, but I believe you'll be better served asking that as a separate question. — J. M.'s missing motivation, May 12 '13 at 01:01
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. — InquisitiveInquirer, May 13 '13 at 08:12
The mathematical difficulties of getting even approximate closed-forms solutions for n > 2 bodies are not related to what happens in a numerical calculation. The calculations bypass any consideration of theory. There can be numerical accuracy problems when two bodies get very close to each other, but in fact, numerical solutions can provide very high accuracy if done with the proper attention to details of the calculations. Navigation for satellites and spacecraft are, I believe, all performed numerically. — Ralph Dratman, May 15 '13 at 05:56

score 24 · Accepted Answer · answered May 11 '13 at 15:48

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.

(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.'s missing motivation, May 11 '13 at 15:56
Wouldn't a symplectic integrator approach only be necessary if you're looking to model complex systems? — InquisitiveInquirer, May 11 '13 at 18:53
@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.'s missing motivation, May 16 '13 at 11:06

score 10 · Answer 2 · edited Oct 05 '18 at 19:30

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

Thanks belisarius, I'll have a look at compacting the code. – InquisitiveInquirer May 11 '13 at 14:31 — InquisitiveInquirer, May 11 '13 at 14:31

3D orbits and inaccuracy over time

EDIT:

2 Answers2

Linked