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:

And this is what they look like some time later:

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?

Answer 1

Some frames from my version of the animation:

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.

Answer 2

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