Locating periodic orbits with Mathematica

Let's use the simple effective potential of the classical restricted three-body problem

pot = (1 - μ)/Sqrt[(x[t] + μ)^2 + y[t]^2] + μ/Sqrt[(x[t] + μ - 1)^2 +
y[t]^2] + (x[t]^2 + y[t]^2)/2 + μ*(1 - μ)/2;
H = 2*pot - (ux[t]^2 + uy[t]^2);

μ = 0.108511220580186;
C0 = 3.35;


DifferentialEquations[H_, x00_, y0_, ux0_, uy0_] :=
Module[{Deq1, Deq2, Deq3, Deq4},
Deq1 = x'[t] == ux[t];
Deq2 = y'[t] == uy[t];
Deq3 = ux'[t] == D[pot, x[t]] + 2*uy[t];
Deq4 = uy'[t] == D[pot, y[t]] - 2*ux[t];

{Deq1, Deq2, Deq3, Deq4, x[0] == x00, y[0] == y0, ux[0] == ux0, uy[0] == uy0}]


The initial conditions of an orbit are:

x00 = 1.3; y0 = 0; ux0 = 0;
tmin = 0; tmax = 4;
pot0 = pot /. {x[t] -> x00, y[t] -> y0};
uy0 = -Sqrt[2*pot0 - C0 - ux0^2]


x00 is just a good initial guess. The code should correct this value (and of course also the initial value of uy0 through the energy integral) until we "hit" a periodic orbit.

Now let's use the approach suggested here

Clear[uy0];
fuy0[x0_] := Solve[(H /. {x[t] -> x0, y[t] -> y0, ux[t] -> ux0, uy[t] -> uy0}) == C0, uy0][[1, 1, 2]];

f[xp_, tp_] :=
Module[{xx = x[xp, fuy0[xp]] /. solp, yy = y[xp, fuy0[xp]] /. solp,
uxx = ux[xp, fuy0[xp]] /. solp,
uyy = uy[xp, fuy0[xp]] /.
solp}, {Norm[{xx[tp], yy[tp], uxx[tp], uyy[tp]} - {xx[0], yy[0],
uxx[0], uyy[0]}], Norm[xx[tp] - xx[0]]}]

DE = DifferentialEquations[H, x0, y0, ux0, uy0];
solp = ParametricNDSolve[DE, {x, y, ux, uy}, {t, tmin, tmax}, {x0, uy0},
Method -> "Adams", PrecisionGoal -> 13,
AccuracyGoal -> 13];

ans = NumberForm[
Quiet@FindRoot[f[xp, tp], {{xp, x00}, {tp, 1}},
PrecisionGoal -> 13, AccuracyGoal -> 13], 20];
xper = xp /. ans[[1, 1]];
tper = tp /. ans[[1, 2]];

Print["x_per = ", NumberForm[xper, 20]]
Print["t_per = ", NumberForm[tper, 20]]


Using version 9.0 the proposed codes fails to obtain the correct results.

So, my question is: How can we fix this problem?

Another issue: If someone wants to propose another, more efficient method, for locating the initial conditions and the period of the periodic orbits, this would be greatly appreciated.

The initial conditions of the periodic orbit are: x_0 = 1.303900184464743 (p_y0 is obviously derived from the energy integral), while it's period is t_per = 3.81159928041479. These results obtained, almost instantly, from a code written in standard FORTRAN 77.

Below I include a plot showing the path of the periodic orbit on the configuration $(x,y)$ plane. Obviously, the numerical integration was performed using the FORTRAN code.

GENERAL COMMENT: I refuse to believe that a modern program, such as Mathematica, cannot locate a simple periodic orbit, while a 40-year old code written in standard FORTRAN 77 delivers the output in less than a second!!!

• I increased the tmax to 100 and looked at the solution for x when x0=1.3 as you start with. It doesn't converge on a cycle, but instead shows oscillations of increasing magnitude. So I'd start by visually ensuring that there is in fact a periodic orbit to be found in your system. Jan 21 '18 at 19:22
• @ChrisK It does exist, according to old good FORTRAN! The periodic point is at 1.303900184464743, while it's period is only 3.81159928041479 time units. Jan 21 '18 at 19:25
• Can you edit your post showing a periodic orbit obtained with your code? Because when I run it with those values, it doesn't appear to be periodic. Do you expect the periodic orbit to be dynamically stable? Jan 21 '18 at 19:47
• I removed the constant term in the potential and the memory consumption went down. I have no patience to wait for it to finish though. Jan 21 '18 at 19:53
• @ChrisK See my edit. No, the particular periodic orbit is called Lyapunov orbit, and it is always highly unstable. Jan 21 '18 at 19:55

This is not an answer but comments do not allow for graphics. Let us plot the potential.

pot /. \[Mu] -> 0.108511220580186 /. f_[t] :> f /. x -> 0
Plot[%, {y, 0, 10}][![enter image description here][1]][1]


we get

The long distance behavior is that of a spring. Now, you called that potential "classical restricted three-body problem". That sounds to me like Celestial Mechanics. Are you sure that you did not miss a parenthesis?

pot = (1 - μ)/Sqrt[(x[t] + μ)^2 + y[t]^2]
+ μ/Sqrt[(x[t] + μ - 1)^2 + y[t]^2]
+ 1/2*(*the next factor is in the denominator*)(x[t]^2 + y[t]^2)


Again, just because of no pictures in comments.

I meant to write that the factor after the comment should be in the denominator. Your original post and the current one have that factor in the numerator. It seems odd to me. Compare with the potentials in the link you posted. I plotted the most exotic of those potentials and it yields what I would expect of a gravitational potential.

Vb = (G*Mb)/(2*a)*(ArcSinh[(x[t] - a)*(y[t]^2 + c^2)^(-1/2)] -
ArcSinh[(x[t] + a)*(y[t]^2 + c^2)^(-1/2)]);
Vb /. {G -> 1, Mb -> 1, a -> 1, c -> 1, f_[t] :> f} /. y -> x
Plot[%, {x, 0, 10}]


.

• You are right! When I copied the code from the notebook to the site it got wrong! See my edit. Jan 21 '18 at 20:36
• @Vaggelis_Z So why have the equations not been corrected? It is really nuisance to run code, only later to read comments strongly suggesting the code is not correct, the poster knows that, and has not bothered to correct it. What a waste of time. Jan 21 '18 at 21:57
• @DanielLichtblau, The poster did change the original post but to an equivalent expression of his original one. I grant him that. Implied in my comments is that we cannot compare the FORTRAN answer with the MMA answer because they might be solving different potentials. They both yield periodic motion but the periods are not the same. I hope the op has an answer for us when he comes back from church … it is Sunday over here ;-) Jan 21 '18 at 23:38
• 1630 CST (approximately when I posted is past noon everywhere east of the IDL... In any case, I suspect this is an apples-to-oranges comparison in that the potentials are very likely different in the two cases. There is the numerator/denominator thing and also I wonder if the last term is meant to have a sqrt like the two that precede it. Jan 22 '18 at 16:48

Not an answer, but an extended comment. I am not sure there is a periodic orbit in your implementation of this system in Mathematica.

Plot[Evaluate[x[1.303900184464743, 3.81159928041479][t] /. solp],
{t, 0, 3.81159928041479}]


looks far from periodic.

• This may have to do with what integrator is used. If the integrator isn't symplectic it might violate conservation of energy (for example) and then you won't ever see a closed orbit. Jan 21 '18 at 20:15