Hi Moderation Team,

Thanks for the plaudit! But if we're going to make this a featured post, let's go ahead and add an animation.

Performing standard numerical time evolution along the pendulum C.E.S. allows us to compute a "phase space race":

This image shows trajectories for six energies, with ten approximations at each energy level, using the same coloring scheme where precision increases toward blue.

Especially for the outer higher-energy trajectories, the poorly-approximated green curves race ahead of schedule, finishing before their due time. The blue curves linger a while longer around the critical points, and arrive at the finish line in more reasonable time.

Since the period is defined as the time for one complete cycle around phase space, the energy dependence of the period is very obvious in the animation.

Hi Todd,

Same to you... I was reading through your ABC book, cool idea.

I agree that the code is difficult to follow without pictures, but remind you that the pictures are already out there in numerous references. And in the future we are likely to see more references, with pictures, better explanations, etc.

Actually the real difficulty is with the optimization, but this is also the main bragging point. Whereas only 10 orders of magnitude are given in the Wolfram Demonstration "Exact and Approximate Relativistic Corrections to the Orbital Precession of Mercury", in the above code we can easily compute 20 orders of magnitude.

If you want to see some pictures, we might as well make some.

• Add phase angle dependence, and list trajectories through 20 orders of magnitude

  TrajectoryList = (RExp[#] /. R -> RSet /. v[n_] :> Q^n v[n]) & /@ 
  Range[20];
  

• User Coulomb potential expansion coefficients for v[ _ ]

  KeplerPhaseSpace =  MapThread[#1 /. #2 &, {TrajectoryList, 
  vRep[PotentialList[[2]], 2 + #] & /@ Range[20]}];
  

• List plot functions for 10 Conserved Energy Surface (CES) contours in phase space

  PlotFList = Flatten[KeplerPhaseSpace /. {Q -> Cos[\[Phi]], b -> #/15} & /@ Range[10]];
  

• Plot the CES contours. Precision increases from green to blue.

  Show[MapThread[  PolarPlot[#1, {\[Phi], 0, 2 Pi}, PlotStyle -> #2] &, {PlotFList, 
  Flatten[Table[  RGBColor[0, (20 - #)/20, #/20] & /@ Range[20], {10}]]}], 
  PlotRange -> All, ImageSize -> 800]
  

• CES Coordinates at order-20 precision

  f20 = {# Cos[\[Phi]], # Sin[\[Phi]], 1/2 b^2} &@(TrajectoryList[[20]] /. 
  vRep[PotentialList[[2]], 22] /. Q -> Cos[\[Phi]]);
  

• Place potential function into phase space with zero momentum

  Potential = {x, 0, 1/2 + 1/(2 (1 + x)^2) - 1/(1 + x)}
  

• Plot CES over phase space with potential function as a tangent curve

  Show[  ParametricPlot3D[f20, {\[Phi], 0, 2 \[Pi]}, {b, 0, 10/15}, 
  ImageSize -> 600, Axes -> False, Boxed -> False, PlotRange -> All],
  ParametricPlot3D[Potential, {x, -1/2, 3}], ImageSize -> 1000 ]
  

These are just plots related to the "Rexp" function. With some additional development of the "dt" function, we could even animate motion of the oscillator along a CES.

Happy Holidays !