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[0] == ux0, x[0] == x0, uy[0] == uy0, y[0] == y0}

Then the module for the computation and plot of the orbit

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

  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[[1]];
  yy[t_] = y[t] /. sol[[1]];
  uxx[t_] = ux[t] /. sol[[1]];
  uyy[t_] = uy[t] /. sol[[1]];

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

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[[1]];
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

enter image description here

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?


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[[1]]
xmin = -15; xmax = 15;
ymin = -15; ymax = 15;

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

we get

enter image description here

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

  • $\begingroup$ 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? $\endgroup$
    – Michael E2
    Commented Jan 7, 2015 at 11:32
  • $\begingroup$ @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. $\endgroup$
    – Vaggelis_Z
    Commented Jan 7, 2015 at 11:35
  • 1
    $\begingroup$ 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. $\endgroup$
    – Michael E2
    Commented Jan 7, 2015 at 11:39
  • $\begingroup$ I don't know what Adams is as an integration method, but you might want a symplectic integrator such as leapfrog. $\endgroup$
    – chris
    Commented Jan 7, 2015 at 14:54
  • $\begingroup$ @chris You can use any integrator you like. I need to know to compute an accurate value of x0using an iterative correction method. $\endgroup$
    – Vaggelis_Z
    Commented Jan 7, 2015 at 15:02

1 Answer 1


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:


(* {tmax->2.38841,x0->13.6924} *)

Alternatively, you might be able to use a WhenEvent.

  • $\begingroup$ 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. $\endgroup$
    – Vaggelis_Z
    Commented Jan 10, 2015 at 14:58
  • $\begingroup$ I don't see why not. $\endgroup$
    – Chris K
    Commented Jan 10, 2015 at 16:40
  • $\begingroup$ 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. $\endgroup$
    – Vaggelis_Z
    Commented Jan 12, 2015 at 21:16
  • $\begingroup$ 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? $\endgroup$
    – Vaggelis_Z
    Commented Jan 13, 2015 at 7:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.