Solve the reduced three-body problem?

Questions
, , ,
1

Hello,
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. 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)

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

u = 0.5

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]

In[184]:= 
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} ]

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



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

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

In[111]:= 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  *)

In[182]:= Plot[Evaluate[J[t]], {t, 1, 30}]

Programming task 1.nb (1.38 MB)

2

Welcome to Wolfram Community! Please make sure you know the rules: https://wolfr.am/READ-1ST

The rules explain how to format your code properly. If you do not format code, it may become corrupted and useless to other members. Please EDIT your post and make sure code blocks start on a new paragraph and look framed and colored like this.

int = Integrate[1/(x^3 - 1), x];
Map[Framed, int, Infinity]

ScreenShot2017-04-28at9.43.33AM.png
ScreenShot2017-04-28at9.43.33AM.png500Ă—135 18.6 KB

3

Hope this is better.

4

Try PrecisionGoal -> 30, or something like that, in DSolve.

5

Thanks for the idea Gianluca. Surprisingly, this turns my plot of the motion into a straight line that that goes to huge negative values, order 10^7. Do you know of anything else that might be missing? If I keep a high upper bound on t I get the error message:

NDSolve::ndsz: At t == 33.82118505471137`, step size is effectively zero; singularity or stiff system suspected.

How do I deal with a singularity? By increasing the step size maybe?

6

Don’t use floating point constants or initial data such as u=.5, but rational exact constants, as they degrade precision.

7

Ok, thanks. This doesn’t resolve the issue, unfortunately.

8

I have no particular problems with this code:

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];
{xsol, ysol} =
 NDSolveValue[{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}, PrecisionGoal -> 30]
ParametricPlot[{xsol[t], ysol[t]},
 {t, 0, 100}]
j[t_] := 1/2*xsol'[t]^2 + 1/2*ysol'[t]^2 +
   Om[xsol[t], ysol[t], u];
Plot[j[t], {t, 0, 100}]
9

We can add one more option and get a solution on which the invariant J will be conserved with high accuracy

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] == 1/2}, {x, y, x', y'}, {t, 
    100}, WorkingPrecision -> 30, MaxSteps -> 10^6];

J[t_] := ((x'[t]^2 + y'[t]^2)/2 + Om[x[t], y[t], u]) /. sol

{ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 100}, 
  PlotRange -> All, AxesLabel -> {"x", "y"}], 
 Plot[J[t] - J[0], {t, 0, 100}, 
  AxesLabel -> {"t", "\[CapitalDelta]J"}]}

3Body.png
3Body.png962Ă—264 102 KB