Message Boards Message Boards

GROUPS:

Solve the reduced three-body problem?

Posted 1 month ago
369 Views
|
8 Replies
|
2 Total Likes
|

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:
8 Replies

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 1 month ago

Hope this is better.

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

Posted 1 month 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?

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

Posted 1 month ago

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

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

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

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