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.


2 Answers 2


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

enter image description here

enter image description here

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[0] == y[0] == 0,
        x'[0] == y'[0] == 1/2}, {x, y}, {t, 100}, MaxSteps -> 500000, 
      Method -> {"Projection", Method -> "StiffnessSwitching", 
        "Invariants" -> J[t]}]

enter image description here

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$:

enter image description here

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.

  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Sep 7, 2018 at 12:01
  • $\begingroup$ 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. $\endgroup$
    – SPPearce
    Sep 7, 2018 at 12:49
  • $\begingroup$ So it depends on what the use-case of the simulation is as to whether the accuracy of one particular simulation is important. $\endgroup$
    – SPPearce
    Sep 7, 2018 at 12:58
  • $\begingroup$ But, why not increase WorkingPrecision until the answer stops changing? $\endgroup$
    – bbgodfrey
    Sep 7, 2018 at 13:01
  • $\begingroup$ 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. $\endgroup$
    – SPPearce
    Sep 7, 2018 at 13:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.