Advanced Runge-Kutta vs Symplectic integrators

230 Views Asked by At

Symplectic integrators are build from compositions of discrete symplectomorphisms. Symplectic integrators do not conserve energy, but energy error is bounded. Since non-linear time translations are generaly not symplectomorphisms, adaptive step sizes do not usually go well with symplectic methods (attempts exist to mix them and are usually cumbersome, e.g. this). Below I consider conservative Hamiltonian systems exclusively.

All comparisons between Runge-Kutta (RK) and symplectic integrators (SI) I have ever seen are carried out between elementary methods, usually of 2nd or 3rd order and always with fixed time steps for RK. It's nice for a textbook to show leapfrog has closed orbits while RK4 goes crazy for a simple Kepler problem, but seeing comparisons like this even in papers is off-putting. I have never seen comparisons between e.g. a 4th order Blanes and Moan's against Tsitouras5(4) using a proportional-integral controller for adaptive step-sizing. I have done some comparisons myself and found out that the energy drift in PI-controlled RK schemes is very similar to the one in SIs, either oscillating in a bounded way or growing extremely slowly as a function of time. If the tolerances in RK are decreased or a higher-order method is used, energy growth can be kept within the oscillation bounds in an equivalent SI scheme. I then thought:

  1. Many RK methods have embedded dense-output formulas and allow for producing a final continuous solution of the same order as the algorithm used. This is a big advantage for event handling, and I don't think SIs have something equivalent.

  2. Controllers pick very large steps outside of problematic zones where SIs, due to the fixed step sizes, are forced to go as slowly as they go near critical points or discontinuities. This is very inefficient.

  3. Each full step in a high-order SI takes dozens of function evaluations, which is a massive problem for systems with many degrees of freedom where functions are slow. Each step in a 5(4) FSAL RK, on the other hand, requires only 5 function evaluations, and the error estimates come from one single function call.

In my head advanced RK methods would then have an advantage, but then again I started to wonder what could go wrong:

  1. I have no idea of how SIs deal with the problem of stiffness - I don't even know if this problem exists for them. For RK methods one can employ automatic stiffness detection and switch to a low-order implicit or RADAU scheme when needed, which are slower than explicit RK methods and, if stiffness is frequent, I assume a SI can be cheaper... But again, I have no idea of how SIs are impacted by stiffness. If it is a problem for them too, then this is not a true issue.

  2. If solutions are always close to problematic regions, i.e. if the Hamiltonian is a very singular function (think of the solar system or a more complicated many-body problem with Coulombic interactions), then the controller will pick tini-tiny steps the whole time and the advantage of adaptive steps vanishes. Of course, it might be that a SI can only deal with the same problem for a fixed step size which is bounded from above by the mean RK step size, so an SI would also perform poorly. I hope SI performs better than RK here, otherwise people in celestial mechanics have been choosing the wrong integrators since the 17th century.

Those things kept going around in my head and, if it hasn't become clear until now, I lean towards thinking that RK algorithms using a smart controller offer an advantage over SIs. I would highly appreciate if someone can confirm or disprove this line of thought with either personal experience or references.

P.S.: I don't believe this question is a duplicate. I have read most other questions about RK vs SI and they don't really touch any of the points I'm raising here.