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

Question

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[0] == y[0] == 0,
    x'[0] == y'[0] == 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][[1]]

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

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

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

(* 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.

Answer 1

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[0] == y[0] == 0,
    x'[0] == y'[0] == 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.

Answer 2

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[0] == y[0] == 0,
        x'[0] == y'[0] == 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.