# Solve the reduced three-body problem?

Posted 3 months ago
630 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
Sort By:
Posted 3 months ago
 Welcome to Wolfram Community! Please make sure you know the rules: https://wolfr.am/READ-1STThe 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] 
Posted 3 months ago
 Hope this is better.
Posted 3 months ago
 Try PrecisionGoal -> 30, or something like that, in DSolve.
Posted 3 months 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 3 months ago
 Don't use floating point constants or initial data such as u=.5, but rational exact constants, as they degrade precision.
Posted 3 months 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"}]} `