Celestial mechanics and Runge Kutta methods

Question

I am working on an example here to simulate the orbit of Earth for one year. As you can see in the notebook, RK45 doesn't conserve energy, and after one simulated year it has spiraled in substantially. That's as expected.

I also don't expect RK23 to conserve energy; however, I have seen it do better with this kind of system, so I gave it a try.

It turns out to do a lot better: with fewer total function evaluations, it loses much less energy and ends up pretty close to the start.

Does anyone know why RK23 seems to do unreasonably well for this example?

I am using the SciPy function solve_ivp, if you want more detailed info on the implementations: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

[And just to be clear, I know that there are other methods that conserve energy; that's not what this question is about.]

Answer 1

According to the documentation for the SciPy function solve_ivp, RK23 is based on the Bogacki-Shampine method, which is implemented in the MATLAB function ode23. Below are numerical results obtained from applying ode23 to a long-time integration of two Hamiltonian systems: a simple double-well example and the OP's earth orbit example.

Double-Well Example

Here I run ode23 on a double-well Hamiltonian system with Hamiltonian function $H(q,p) = (1/2) p^2 + (1/4) (q^2 - 1)^2$ using the following MATLAB code

ff=@(t,y) [y(2); y(1)-y(1).^3];   
opts=odeset('RelTol',tol); 
[t,y]=ode23(ff,[0 5000],[0; 2],opts);

for the tol values indicated in the figure titles below starting with the default relative tolerance of $10^{-3}$. Basically this code numerically integrates $$ \dot{q} = p \;, \quad \dot{p} = q - q^3 \;, \quad (q(0),p(0)) = (0,2) \;, $$ for a (long) time span of $5000$. ode23 outputs a vector of times $t$ and a matrix $y$ whose rows are the corresponding numerical approximation.

Below are the outputted discrete trajectories in phase space. The dots are made a bit lighter with time. The solid red curve is the level set of the Hamiltonian corresponding to the initial point. The actual solutions lie on this red curve for all time since they preserve $H$. In contrast, the outputted dots seem to converge to the right well, and the ones with smaller tol seem to take longer to converge.

OP's Earth Orbit Example

Following the OP's description, one can similarly simulate the earth revolving around the sun using

m1=1.989e30;
m2=5.972e24;
G=6.674e-11;
ff=@(t,y) [y(3); y(4); ...
  -G*m1*y(1)/(y(1)^2+y(2)^2)^(3/2); ...
  -G*m1*y(2)/(y(1)^2+y(2)^2)^(3/2)];
T=100*31556925.9747; % time-span is 100 years!
opts=odeset('RelTol',tol);
[t,y]=ode23(ff,[0 T],[147*1e9; 0; 0; -30300],opts);

The following figure shows the relative energy error after 100 years for two different tol values.

Discussion

Since we see a systematic energy drift in both of these relatively simple Hamiltonian test problems, these numerical counterexamples illustrate that ode23 is probably not a geometric integrator. For example, a geometric integrator that preserves the symplectic form of the Hamiltonian system (called symplectic integrators), typically do not preserve energy, but they do have bounded energy errors over long-time simulations. One can construct adaptive geometric integrators, but this is a bit tricky. See, e.g.,

Calvo, M. P.; López-Marcos, M. A.; Sanz-Serna, J. M., Variable step implementation of geometric integrators, Appl. Numer. Math. 28, No. 1, 1-16 (1998). ZBL0930.65136..