Try TaylorIntegration.jl for integrators of arbitrarily high order. Cc @lbenet
Thank you, Henon-Heiles is trivial, I have that working with and without symplectic integrators. The circular restricted three-body problem however is nasty:
The usual “synodic” coordinate system is not an inertial system, it is rotating to seperate the motion of the heavy masses off. The heavy masses are at rest making the ODE for the test particle autonomous and the problem then has only four degrees of freedom. This however has a consequence:
Energy is not an integral of motion ! Also, angular momentum is not an integral of motion. The conserved quantity is the Jacobi Integral which is basically the sum of energy and angular momentum. There exists a second integral (conserved quantity) but that cannot be expressed analytically as far as I am aware. Further, if you look at the formulation as a Newtonian equation of motion for the test particle (so second order differential equations) you see that the force is velocity dependent. As a consequence common symplectic integrators, e.g Verlets or Forrest-Ruth type, make unfortunately no sense as they conserve the wrong quantity and are not applicable anyway due to the velocity dependence of the force 
I originally was planning this for a lecture demonstration (they did that on 1960s computers, it must be simple, right ?) and in fact it is a nice problem where many things fail. However, without very good preparation and introduction … I will probably pass it and what remains is my personal curiosity.
I am already using a callback to monitor the Jacobi integral and abort if it deviates too much. Eventually I can play with this. My observation is however that small deviations of the integral (1e-12) already for the unstable orbits completelely change the orbit.
Remark: Of course one could work in cartesian coordinates and an inertial system and integrate alle masses (so 18 degrees of freedom). However, near collision orbits still need regularization as far as I read on this.
Yes, Eduard Stiefels work. I know that. I gave up early. Probably everything is in there if one can find and specialise it from the general case
If you are still interested this has a executive level summary. Discussion of regularization is missing, though.
That method beyond Euler is new to me as a practical one. The general idea is of course simple enough. From a first look there might be a superficial similarity to the recursion approach. However, if floating point rounding error is the actual problem this will not help.
Addendum:
Reading this article
and looking at the Van der Pol example I am quite sure that NASAs
recursion method is a very close relative, see “Power series solutions
of the Thiele-Burrau regularized planar restricted three-body
problem”. The article
also contains an application to the non-planar RTBP. Sections 5.1.3 and
5.1.4 discuss floating point arithemtics effects. Quote “Both graphics
show that the error behaviour seem to be dominated by the noise of the
floating point arithmetic.”, so that appears to be an issue ?
I’d just like to mention that there’s an example in the TaylorIntegration.jl docs where the planar RTBP is integrated. Our intention in that example was to show the interoperability between TaylorIntegration.jl and the DifferentialEquations.jl ecosystem. The example shows that TaylorIntegration.jl is able to preserve the Jacobi integral to \sim 10^{-14} (relative), despite having several close encounters with the secondary. Maybe that example can be of some help.
Of course, when dealing with really close encounters (e.g., if the test particle enters the Hill radius of either primary), then regularization helps a lot to preserve the integrals of motion. Some years ago, I used A. Jorba’s Taylor package (in C) to integrate the planar RTBP and was able to get nice accuracies (i.e., close to the machine epsilon) when dealing with near-collision orbits with the secondary, with Levi-Civita regularization. What I interpret from the quote you posted from Jorba’s paper, is that Taylor’s method is able to reach, at each step, the roundoff accuracy of 64-bit floating point, which is near the best you can ask for with such arithmetic.
Now, when dealing with BigFloats in Julia, with TaylorIntegration.jl you can adjust the parameters order and abstol in taylorinteg, so that you are able to reach the numerical accuracy of BigFloat. There’s a simple (yet, not very user-friendly) test in TaylorIntegration.jl where we actually check this. There, we compute the complete elliptic integral of the 1st kind K(k) to an accuracy of about \sim 10^{-77} using BigFloats, which we then use to compute the solution to the simple pendulum to the same accuracy.
It would be interesting to implement the Thiele-Burrau regularization for the RTBP with TaylorIntegration.jl and BigFloats and see how it performs for your problem; if there’s still some interest I’d be happy to give it a try!
1 Like