# How to locate the position of a periodic orbit

These are the equations of the dynamical system

Vn = (-G*Mn)/Sqrt[x[t]^2 + y[t]^2 + cn^2];
Vd = (-G*Md)/Sqrt[x[t]^2 + y[t]^2 + (s + h)^2];
Vh = (-G*Mh)/Sqrt[x[t]^2 + y[t]^2 + ch^2];
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)]);
pot = Vn + Vd + Vh + Vb;
H = 1/2*(ux[t]^2 + uy[t]^2) + pot - om*(x[t]*uy[t] - y[t]*ux[t]);


and these are the values of the parameters

G = 1; Mn = 400; cn = 0.25;
Md = 7000; s = 3; h = 0.175;
Mb = 3500; a = 10; c = 1;
Mh = 20000; ch = 20;
om = 4.5;
H0 = -3180;


The initial conditions of the orbit are

x00 = 10.77; y0 = 0; ux0 = 0;
Ht = H /. {x[t] -> x00, y[t] -> y0, ux[t] -> ux0};
pot0 = pot /. {x[t] -> x00, y[t] -> y0};
py0 = x00*om - Sqrt[x00^2*om^2 + 2*(H0 - pot0)];
sol = Solve[Ht == H0];
uy0 = uy[t] /. sol[]
tmin = 0; tmax = 1;


The set of the equations of motion

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

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


and the numerical integration

DE = DifferentialEquations[H, om, x00, y0, ux0, uy0];
sol = NDSolve[DE, {x[t], y[t], ux[t], uy[t]}, {t, tmin, tmax},
MaxSteps -> Infinity, Method -> "Adams",
PrecisionGoal -> 12, AccuracyGoal -> 12];
xx[t_] = x[t] /. sol[];
yy[t_] = y[t] /. sol[];
uxx[t_] = ux[t] /. sol[];
uyy[t_] = uy[t] /. sol[];


For x00 = 10.77 the corresponding orbit is the follwoing

plot = ParametricPlot[{xx[t], yy[t]}, {t, tmin, tmax}, Axes -> False,
Frame -> True, AspectRatio -> 1, PlotStyle -> Black,
AspectRatio -> 1, PlotRange -> All] We see that the orbit is not periodic. However if we use x00 = 10.77403 we get which is indeed a periodic orbit.

My question is obviously the following: how can I locate the exact (let's say with 10 decimal digits) position of the periodic orbit? Somehow inside the NDSolve there should be an iterative process changing the value of x00 until it hits the periodic point.

The corresponding FORTRAN code indicates that the position of the periodic orbit is at x00 = 10.774029735833850. So any provided method here must give the same result.

NOTE: The energy level H0 = -3180 should be remain the same while searching for the x00 value of the periodic orbit. x00 is always in the interval [9,12], so the initial guess 10.77 should be corrected somehow so as to hit the exact the periodic point. Also for x = x00 it should be y0 = ux0 = 0.

EDIT

DO loop for variable value of the energy

data = {};
Do[
x00 = 10.5; y0 = 0; ux0 = 0;
tmin = 0; tmax = 1;
Ht = H /. {x[t] -> x00, y[t] -> y0, ux[t] -> ux0};
pot0 = pot /. {x[t] -> x00, y[t] -> y0};
py0 = x00*om - Sqrt[x00^2*om^2 + 2*(H0 - pot0)];
sol = Solve[Ht == H0];
uy0 = uy[t] /. sol[];
Clear[uy0];
fuy0[x0_] :=
Solve[(H /. {x[t] -> x0, y[t] -> y0, ux[t] -> ux0, uy[t] -> uy0}) ==
H0, 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,
yy, uxx, uyy}], Norm[xx[tp] - xx]}]
DE = DifferentialEquations[H, om, x0, y0, ux0, uy0];
solp = ParametricNDSolve[
DE, {x, y, ux, uy}, {t, tmin, tmax}, {x0, uy0},
MaxSteps -> Infinity, Method -> "Adams", PrecisionGoal -> 12,
AccuracyGoal -> 12];
pos = Quiet@
FindRoot[f[xp, tp], {{xp, x00}, {tp, .5}}, PrecisionGoal -> 12,
AccuracyGoal -> 12];
xper = xp /. pos[];
tper = tp /. pos[];
AppendTo[data, {xper, tper}],
{H0, -3180, -3170, 1}
]

• The answer to The Orbit and Perigee of the Flamsteed comet should be helpful. – Artes Oct 6 '15 at 10:42
• Are all of {x0, y0, ux0, uy0} free parms? – Dr. belisarius Oct 6 '15 at 13:03
• @belisarius x0 is to be determined, y0 = 0, ux0 = 0, uy0 ---> obtained from the energy integral. BTW nice to hear from you again! – Vaggelis_Z Oct 6 '15 at 13:08
• One thing to point out, numerical solvers often exhibit energy drift where the system appears to lose or gain energy as the simulation progresses. So, any fixed point you find needs to be considered suspicious until you've demonstrated that the simulation does not exhibit drift. – rcollyer Oct 6 '15 at 13:18
• @rcollyer Yes, any numerical integrator has an error but in this case the error is beyond the twelfth significant figure. – Vaggelis_Z Oct 6 '15 at 13:20

With uy0 defined in terms of x0 as

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


the criterion for a repeated orbit as

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, yy, uxx, uyy}],
Norm[xx[tp] - xx]}]


and other quantities as in the question, then

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

NumberForm[Quiet@
FindRoot[f[xp, tp], {{xp, x00}, {tp, .5}}, PrecisionGoal -> 12, AccuracyGoal -> 12],
15]
(* {xp -> 10.774029731533837, tp -> 0.5320581303031949} *)


where the first number is the x0 initial condition, and the second number the period. The calculation is virtually instantaneous.

Clear[xx, yy, uxx, uyy];
xx = x[xp, fuy0[xp]] /. ans[] /. solp;
yy = y[xp, fuy0[xp]] /. ans[] /. solp;
uxx = ux[xp, fuy0[xp]] /. ans[] /. solp;
uyy = uy[xp, fuy0[xp]] /. ans[] /. solp;
plot = ParametricPlot[{xx[t], yy[t]}, {t, tmin, tmax}, Axes -> False,
Frame -> True, AspectRatio -> 1, PlotStyle -> Black, AspectRatio -> 1, PlotRange -> All] Response to Edit with new code

Vn = (-G*Mn)/Sqrt[x[t]^2 + y[t]^2 + cn^2];
Vd = (-G*Md)/Sqrt[x[t]^2 + y[t]^2 + (s + h)^2];
Vh = (-G*Mh)/Sqrt[x[t]^2 + y[t]^2 + ch^2];
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)]);
pot = Vn + Vd + Vh + Vb;
H = 1/2*(ux[t]^2 + uy[t]^2) + pot - om*(x[t]*uy[t] - y[t]*ux[t]);
G = 1; Mn = 400; cn = 0.25; Md = 7000; s = 3; h = 0.175; Mb = 3500; a = 10;
c = 1; Mh = 20000; ch = 20; om = 4.5;
x00 = 10.77; y0 = 0; ux0 = 0; tmin = 0; tmax = 1;
DifferentialEquations[H_, om_, x00_, y0_, ux0_, uy0_] :=
Module[{Deq1, Deq2, Deq3, Deq4},
Deq1 = x'[t] == ux[t] + om*y[t];
Deq2 = y'[t] == uy[t] - om*x[t];
Deq3 = ux'[t] == -D[pot, x[t]] + om*uy[t];
Deq4 = uy'[t] == -D[pot, y[t]] - om*ux[t];
{Deq1, Deq2, Deq3, Deq4, x == x00, y == y0, ux == ux0, uy == uy0}];
data = {};
Do[Clear[uy0];
fuy0[x0_] := Solve[(H /. {x[t] -> x0, y[t] -> y0, ux[t] -> ux0, uy[t] -> uy0}) ==
H0, uy0][[1, 1, 2]];
DE = DifferentialEquations[H, om, x0, y0, ux0, uy0];
solp = ParametricNDSolve[DE, {x, y, ux, uy}, {t, tmin, tmax}, {x0, uy0},
MaxSteps -> Infinity, Method -> "Adams", PrecisionGoal -> 12, AccuracyGoal -> 12];
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, yy, uxx, uyy}],
Norm[xx[tp] - xx]}];
ans = NumberForm[Quiet@FindRoot[f[xp, tp], {{xp, x00}, {tp, .5}},
PrecisionGoal -> 12, AccuracyGoal -> 12], 15];
xper = xp /. ans[[1, 1]];
tper = tp /. ans[[1, 2]];
AppendTo[data, {xper, tper}], {H0, -3180, -3170, 1}]

data
(* {{10.774, 0.532058}, {10.7705, 0.53089}, {10.7668, 0.529734}, {10.7631, 0.52859},
{10.7594, 0.527458}, {10.7556, 0.52634}, {10.7517, 0.525235}, {10.7478, 0.524144},
{10.7439, 0.523068}, {10.7399, 0.522006}, {10.7358, 0.52096}} *)

• You are a life saver! I owe you a six-pack! – Vaggelis_Z Oct 6 '15 at 15:49
• Just a comment. If I want to directly plot the periodic orbit after calculating the position what should I do? – Vaggelis_Z Oct 6 '15 at 15:52
• @Vaggelis_Z I have added the plot and associated code, as you requested. Thanks for accepting the answer. Best wishes. – bbgodfrey Oct 6 '15 at 16:07
• Please see my EDIT. I created a loop for finding the position of the orbits for variable value of the energy. However it does not work. – Vaggelis_Z Oct 6 '15 at 16:15
• @Vaggelis_Z Remember that FindRoot requires reasonable guesses. For H0 = -2730, choosing x00 = 2.4 yields {xp -> 2.4343843379119816, tp -> 0.3677217902335441}. Eventually, you will have to change 0.5 in FindRoot too. Also, remember that there may be more than one solution for a given H0. – bbgodfrey Oct 6 '15 at 20:35

## Edit

Here is a solution that addresses energy level:

(* parameter dep. system *)
DE = DifferentialEquations[H, om, X, Y, UX, UY] ;
(* function of initial coordinates for fixed end time *)
f1[val_?NumberQ] := With[
{T=val},
ParametricNDSolveValue[
DE,
{x[T],y[T],ux[T],uy[T]},
{t, 0, T},
{X,Y,UX,UY},
MaxSteps -> Infinity,
Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 20, "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients", "ImplicitSolver" -> {"Newton", AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision, "IterationSafetyFactor" -> 1}},
WorkingPrecision->MachinePrecision
]] ;
(* find fixed point for fixed end time *)
f2[val_?NumberQ] := With[
{fun = f1[val]},
{X,Y,UX,UY} /. FindRoot[fun[X,Y,UX,UY]=={X,Y,UX,UY},{X,x00},{Y,y0},{UX,ux0},{UY,uy0},Evaluated->False]
] ;
(* value of Hamiltonian for given fixed point *)
f4[t_?NumberQ] := Apply[f3,f2[t]] ;
(* find period *)
per = Q /. FindRoot[H0 - f4[Q] == 0,{Q,1}]
(* recover f.p. *)
fp = f2[per]
(* check Hamiltonian *)


1.06412 (* period *)

{10.7739, 0.0210223, 0.0700877, 37.2} (* initial condition *)


You need to solve a fixed point problem $\varphi(x) = x$ where $x$ is a vector of initial values and $\varphi$ is a solution at $t=1$.

First, define a parameter dependent system:

DE = DifferentialEquations[H, om, X, Y, UX, UY] ;
f = ParametricNDSolveValue[
DE,
{x[tmax],y[tmax],ux[tmax],uy[tmax]},
{t, tmin, tmax},
{X,Y,UX,UY},
MaxSteps -> Infinity, Method -> "Adams", PrecisionGoal -> 12, AccuracyGoal -> 12
] ;


Then solve f.p. problem:

{xi,yi,uxi,uyi} = {X,Y,UX,UY} /. FindRoot[f[X,Y,UX,UY]=={X,Y,UX,UY},{X,x00},{Y,y0},{UX,ux0},{UY,uy0},Evaluated->False]


DE = DifferentialEquations[H, om, xi, yi, uxi, uyi] ;
sol = NDSolve[DE, {x[t], y[t], ux[t], uy[t]}, {t, tmin, tmax},
MaxSteps -> Infinity, Method -> "Adams",
PrecisionGoal -> 12, AccuracyGoal -> 12];
xx[t_] = x[t] /. sol[];
yy[t_] = y[t] /. sol[];
uxx[t_] = ux[t] /. sol[];
uyy[t_] = uy[t] /. sol[];
{xi,yi,uxi,uyi}
{x[t],y[t],ux[t],uy[t]} /. sol /. t -> tmax // Flatten
plot = ParametricPlot[{xx[t], yy[t]}, {t, tmin, tmax}, Axes -> False,
Frame -> True, AspectRatio -> 1, PlotStyle -> Black,
AspectRatio -> 1, PlotRange -> All]

• Something is wrong here. Your solution reports that the periodic point is at x00 = 10.552099826486845 which is wrong. – Vaggelis_Z Oct 6 '15 at 8:45
• It is very confusing that you talk about periodic orbit while referring only to $x$. The solution given provides values for all phase space coordinates. Have you checked the value of $x$ at $t=2$ for your (FORTRAN) initial condition? Does it stay fixed? – I.M. Oct 6 '15 at 9:44
• The periodic orbit is unstable, so the orbit is no longer periodic after a couple of loops. I want to locate the x00 position after a time interval of one period which means that the orbit has performed only one loop. Does your solution give x00 = 774029735833850 after t = 0.5320581245422367? – Vaggelis_Z Oct 6 '15 at 9:50
• Is the time at which you expect the orbit to overlap with its initial condition fixed by the system? If I read correctly, I.M.'s solution looks for the case where the state at time t=1 is equal to the initial state. But your example completes its loop at around t=0.53 – Jason B. Oct 6 '15 at 9:53
• I believe event location would be a less computationally intensive strategy here... why not try with WhenEvent[]? – J. M. will be back soon Oct 6 '15 at 13:38