# 3D orbits and inaccuracy over time

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?

-
How about adding PerformanceGoal -> "Quality"... – cormullion May 11 '13 at 12:40
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. en.wikipedia.org/wiki/N-body_problem – user7388 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
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 5:56

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[
c {x[t], y[t], z[t]}/Norm[{x[t], y[t], z[t]}]^3],
{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. 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
@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}]
`
-
Thanks belisarius, I'll have a look at compacting the code. – user7388 May 11 '13 at 14:31