# Locating Periodic Orbits

Here is the code for the numerical integration of an orbit. First the module for the definition of the equations of motion.

DifferentialEquations[H_, om_, x0_, y0_, ux0_, uy0_] :=
Module[{Deq1, Deq2, Deq3, Deq4},
Deq1 = x'[t] == D[H, ux[t]];
Deq2 = y'[t] == D[H, uy[t]];
Deq3 = ux'[t] == -D[H, x[t]];
Deq4 = uy'[t] == -D[H, y[t]];

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


Then the module for the computation and plot of the orbit

PlotTrajectory[H_, om_, x0_, y0_, ux0_, uy0_, tmin_, tmax_, xmax_] :=
Module[{},

DE = DifferentialEquations[H, om, x0, 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[];

plot = ParametricPlot[{xx[t], yy[t]}, {t, tmin, tmax},
Axes -> False, Frame -> True, AspectRatio -> 1,
PlotStyle -> Black, AspectRatio -> 1, PlotRange -> xmax];
Show[plot]
]


The definition of the Hamiltonian

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];
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 + Vb;
H = 1/2*(ux[t]^2 + uy[t]^2) + pot - om*(x[t]*uy[t] - y[t]*ux[t]);


The values of the involved parameters

G = 1; Mn = 400; cn = 0.25;
Md = 7000; s = 3; h = 0.175;
Mb = 3500; a = 4; c = 1;
om = 1.25;
H0 = -900;


The initial conditions of the orbit

x0 = 13.7087; y0 = 0; ux0 = 0;
tmin = 0; tmax = 50;
Ht = H /. {x[t] -> x0, y[t] -> y0, ux[t] -> ux0};
sol = Solve[Ht == H0];
uy0 = uy[t] /. sol[];
xmin = -15; xmax = 15;
ymin = -15; ymax = 15;


and finally the execution

S1 = PlotTrajectory[H, om, x0, y0, ux0, uy0, tmin, tmax, xmax]


which yields We see that the orbit is almost periodic. What I want is to compute the exact position of the periodic orbit. When I say exact I mean accurate enough with 10 decimal digits precision. This orbit has y0 = ux0 = 0, so only x0 should be computed. The value 13.7087 is just an initial good guess. I want to find a way to scan x-axis around this initial value in an appropriate interval and calculate p_y from the value of the total energy. Let the trajectory run for a complete revolution and check whether it returns to the initial point. If not, then change x0 a little to improve the point of return and repeat until the orbit returns to the initial point within the numerical precision we want. So, what is needed is an iterative process for correction of x0 until the periodic point is found.

Any suggestions?

EDIT

For a = 3 Chris K's method gives initial condition x0 = 14.37555820137562, However if I integrate the orbit

x0 = 14.37555820137562; y0 = 0; ux0 = 0;
tmin = 0; tmax = 50;
Ht = H /. {x[t] -> x0, y[t] -> y0, ux[t] -> ux0};
sol = Solve[Ht == H0];
uy0 = uy[t] /. sol[]
xmin = -15; xmax = 15;
ymin = -15; ymax = 15;

S1 = PlotTrajectory[H, om, x0, y0, ux0, uy0, tmin, tmax, xmax]


we get which is far from periodic orbit. So, what is going on?

• On the one hand you state that you want the exact position and on the other that you want an iterative process for approximating it. Do you mean the latter? Jan 7 '15 at 11:32
• @MichaelE2 The iterative process is needed so as to obtain the exact value of x0 by correcting little by little the initial guess. When I mean exact I mean of course a good approximation let's say 10 decimal digits. Jan 7 '15 at 11:35
• Since "exact" tends to be interpreted with a specific meaning in Mathematica that is different than your intention, perhaps "accurate" or "accurate to a user-specified error, such as 10 digits" might be a clearer way to describe the position you seek. Jan 7 '15 at 11:39
• I don't know what Adams is as an integration method, but you might want a symplectic integrator such as leapfrog. Jan 7 '15 at 14:54
• @chris You can use any integrator you like. I need to know to compute an accurate value of x0using an iterative correction method. Jan 7 '15 at 15:02

You could use FindRoot to tweak the period (tmax) and x0 to get back to the starting point, with the other coordinates of the system fixed:

map[H_,om_,x0_?NumericQ,y0_,ux0_,uy0_?NumericQ,tmin_,tmax_?NumericQ,xmax_]:=Module[{},
DE=DifferentialEquations[H,om,x0,y0,ux0,uy0];
sol=NDSolve[DE,{x,y,ux,uy},{t,tmin,tmax},MaxSteps->Infinity,
Return[Evaluate[{x[tmax]/.sol,uy[tmax]/.sol}]]
]

Clear[x0,tmax];
FindRoot[map[H,om,x0,y0,ux0,uy0,tmin,tmax,xmax]=={x0,uy0},{tmax,2.3},{x0,13.7038}]
(* {tmax->2.38841,x0->13.6924} *)


Alternatively, you might be able to use a WhenEvent.

• The numerical method you proposed works perfectly. Now is it possible to create a DO loop for computing the periodic point for variable a let's say when a is between 1 and 9. In other words create a list containing three columns a, x0, tper. Jan 10 '15 at 14:58
• I don't see why not. Jan 10 '15 at 16:40
• Please see my EDIT. I cannot figure out how exactly your module map works but it seems that something is wrong. Does your module keep the energy fixed at -900? Anyway, why for a = 3 the reported x0 is not the periodic point? Please be so king and take a look at this. Jan 12 '15 at 21:16
• UPDATE: Your module works but it does not hold the energy constant. I took the computed x00 and uy0 and inserted them into the Hamiltonian like H00 = H /. {x[t] -> x00, y[t] -> 0, ux[t] -> 0, uy[t] -> uy0}. The value is H00 = -873.553 not -900 as it should be. So yes the code computes a periodic orbit but not for the desired energy level. How can this be fixed? Jan 13 '15 at 7:55