# Proper use of NDSolve in the context of the reduced 3 body problem

I am currently trying to solve the reduced three body problem in Mathematica. I have the equations of motion

$\ddot{x}-2\dot{y}=-\frac{\partial\Omega}{\partial x }$

$\ddot{y}+2\dot{x}=-\frac{\partial\Omega}{\partial y }$

$\Omega=-\frac{1}{2}\mu r_1^2-\frac{1}{2}(1-\mu)r_2^2-\frac{\mu}{r_1}-\frac{1-\mu}{r_2}$

Where $r_{1,2}$ are the distances between the small body and the two massive bodies , and $\mu$ is such that the ratio of the masses is $\mu:1-\mu$. I have shown that $J=\frac{1}{2}\dot{x}^2+\frac{1}{2}\dot{y}^2+\Omega(x,y)$ is conserved, by differentiating wrt time. Then to obtain the motion I wrote the code below, but I cannot seem to get $J$ to be conserved. Does anyone have a clue why?

Plotting $x(t)$ and $y(t)$ gives some discontinuous derivatives at $t=6s$ which is when there is the first major change in the value of $J$. I tried forcing NDSolve to use Runge-Kutta but that didn't help, forcing it to use Euler gave nonsense (outward spiralling circular orbit, as if there were no gravity)

r1[x_, y_, u_] := Sqrt[ (x + 1 - u)^2 + y^2]

u = 1/2

r2[x_, y_, u_] := Sqrt[(x - u)^2 + y^2]

Om[x_, y_, u_] := -1/2*u*r1[x, y, u]^2 - 1/2*(1 - u)*r2[x, y, u]^2 -
u/r1[x, y, u] - (1 - u)/r2[x, y, u]

sol = NDSolve[ {x''[t] - 2 y'[t] == -D[Om[x[t], y[t], u], x[t]],
y''[t] + 2 x'[t] == -D[Om[x[t], y[t], u], y[t]], x == y == 0,
x' == y' == 0.5}, {x, y}, {t, 100} ]

ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 100},
PlotRange -> All]
J[t_] := 1/2*dX[t]^2 + 1/2*dY[t]^2 + Om[X[t], Y[t], u]

pX[t_] := x[t] /. sol
X[t_] := pX[t][]

pdX[t_] := x'[t] /. sol
dX[t_] := pdX[t][]

pY[t_] := y[t] /. sol
Y[t_] := pY[t][]

pdY[t_] := y'[t] /. sol
dY[t_] := pdY[t][]

(* dX[t] gives derivative of x[t], X[t] gives x[t]. We introduce X
and Y to obtain reals, not one element lists, which correspond, e.g.,
to x[t]/.sol  *)

Plot[Evaluate[J[t]], {t, 1, 30}]


I have also tried modifying PrecisionGoal in NDSolve but this leads to straight line trajectories (no gravity again) for PrecisionGoal greater than about 20.

One way I am trying to verify the accuracy of the simulation is by generating the Moulton orbits described in this paper: http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?db_key=AST&bibcode=1967AJ.....72..373S. This has been achieved.

The issue is, in fact, WorkingPrecison. Whenever r1 or r2 approach zero, high accuracy is needed for the computation. Try

sol = NDSolve[ {x''[t] - 2 y'[t] == -D[Om[x[t], y[t], u], x[t]],
y''[t] + 2 x'[t] == -D[Om[x[t], y[t], u], y[t]], x == y == 0,
x' == y' == 1/2}, {x, y, x', y'}, {t, 100}, WorkingPrecision -> 75,
MaxSteps -> 500000, Method -> "StiffnessSwitching"];
ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 100}, PlotRange -> All,
ImageSize -> Large, AxesLabel -> {x, y}, LabelStyle -> {Black, Bold, Medium}]
J[t_] = ((x'[t]^2 + y'[t]^2)/2 + Om[x[t], y[t], u]);
Plot[J[t]/.sol, {t, 1, 100}, PlotRange -> All, ImageSize -> Large, AxesLabel -> {t, J},
LabelStyle -> {Black, Bold, Medium}]  WorkingPrecision -> 75 seems to be large enough; 30 definitely is not.

Another thing to try is the Projection method of NDSolve:

u = 1/2;
r1[x_, y_, u_] = Sqrt[(x + 1 - u)^2 + y^2];
r2[x_, y_, u_] = Sqrt[(x - u)^2 + y^2];
Om[x_, y_, u_] = -1/2*u*r1[x, y, u]^2 - 1/2*(1 - u)*r2[x, y, u]^2 -
u/r1[x, y, u] - (1 - u)/r2[x, y, u];
J[t_] = 1/2*x'[t]^2 + 1/2*y'[t]^2 + Om[x[t], y[t], u];

sol2 = NDSolve[{x''[t] - 2 y'[t] == -D[Om[x[t], y[t], u], x[t]],
y''[t] + 2 x'[t] == -D[Om[x[t], y[t], u], y[t]], x == y == 0,
x' == y' == 1/2}, {x, y}, {t, 100}, MaxSteps -> 500000,
Method -> {"Projection", Method -> "StiffnessSwitching",
"Invariants" -> J[t]}] This gives a result in under a second compared to a minute for bbgodfrey's solution, but still has rather large oscillations in the invariant of the interpolating functions. InvariantErrorPlot shows these as smaller, but I had to play with the syntax to get that to work (define extra variables xd,yd in the NDSolve).

Following the discussions in the comment section, note that this system is chaotic, and this method begins to visually diverge with the solution with WorkingPrecision -> 75 at around $t=20$: This can be fixed by increasing the WorkingPrecision with the Projection method, but that takes more computational time. It doesn't require as high WP to get the same results to $t=100$ though, WP=40 instead of WP=75, which takes about a third of the computational time.

Note that changing the initial conditions in the WP=75 solution by $10^{-15}$, the solutions diverge around $t=90$, due to the chaotic nature of the system.

• Nice approach. However, I noticed that your plot differed substantially from mine, so I investigated. With WorkingPrecision -> 75 your plot becomes essentially the same as mine. The point, I believe, is that the "Projection" method does conserve the invariant well, but doing so does not necessarily assure an accurate overall solution. Sep 7 '18 at 12:01
• Yes, I noticed that they were different. They begin to diverge at about $t=20$, and once they do diverge the chaotic nature of the system means they are very different. Sep 7 '18 at 12:49
• So it depends on what the use-case of the simulation is as to whether the accuracy of one particular simulation is important. Sep 7 '18 at 12:58
• But, why not increase WorkingPrecision  until the answer stops changing? Sep 7 '18 at 13:01
• Because it takes 100 times longer. If you change the initial conditions by $10^-15$ with your solution with WorkingPrecision->75, the solutions diverge at $t=90$. If the goal is determining a spacecraft launch then yes, you want as accurate as possible. If the goal is just to plot a plausible orbit then the fast answer is good. Sep 7 '18 at 13:24