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:
