Message Boards Message Boards

Solve the reduced three-body problem?

Posted 6 years ago

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}]
Attachments:
POSTED BY: thomas rialan
8 Replies

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

fig1

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}]
POSTED BY: Gianluca Gorni
Posted 6 years ago

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

POSTED BY: thomas rialan

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

POSTED BY: Gianluca Gorni
Posted 6 years ago

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?

POSTED BY: thomas rialan

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

POSTED BY: Gianluca Gorni
Posted 6 years ago

Hope this is better.

POSTED BY: thomas rialan

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]

enter image description here

POSTED BY: Moderation Team
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract