# 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:
Answer
8 Replies
Sort By:
Posted 1 month 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] 
Answer
Posted 1 month ago
 Hope this is better.
Answer
Posted 1 month ago
 Try PrecisionGoal -> 30, or something like that, in DSolve.
Answer
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?
Answer
Posted 1 month ago
 Don't use floating point constants or initial data such as u=.5, but rational exact constants, as they degrade precision.
Answer
Posted 1 month ago
 Ok, thanks. This doesn't resolve the issue, unfortunately.
Answer
Posted 1 month ago
 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}] 
Answer
Posted 1 month ago
 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"}]} `
Answer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments