High Order Symplectic Integrators (iPython)

Tamayo & Brown (2019) for an overview of these integrators.

Let us start by importing REBOUND, as well as numpy and matplotlib.

import rebound
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt

High order symplectic integrators with a fixed timestep are well suited for planetary systems in which planets orbit the primary mass on almost Keplerian orbits. The planet-planet interactions need to be a perturbation. If they are not, a different integrator, such as IAS15 or MERCURIUS is better suited.

On system that we know is stable is the outer Solar System. So we will use this as our test case. We can either import accurate data from NASA Horizons, or, because this is just a test, use some of the initial conditions which come with REBOUND to setup a simulation.

from rebound import data
sim = rebound.Simulation()
data.add_outer_solar_system(sim) # either this, or add the planets manually
rebound.OrbitPlot(sim);
../_images/HighOrderSymplectic_3_0.png

We’ll integrate the outer Solar System for 1000 years into the future and measure the energy error along the way. We do that at random intervals to make sure we don’t have any aliasing with an orbital period. The following function runs the simulation, and then returns the error.

def measure_energy(sim):
    Nsamples = 1000
    tmax = 2.*np.pi*1e3 # 1000 years
    t_samples = tmax*np.sort(np.random.random(Nsamples))
    E0 = sim.calculate_energy() # initial energy
    Emax = 0. # maximum energy error
    for t in t_samples:
        # we do not want to change the timestep to reach t exactly, thus
        # we need to set exact_finish_time=False and slighlty overshoot.
        sim.integrate(t,exact_finish_time=False)
        E = sim.calculate_energy()
        Emax = max(Emax, np.abs((E-E0)/E0))
    return Emax

To get an idea how our integrators are behaving, we want to run this simulation for various timesteps. So let us set up an array of timesteps from 0.001 to 1 orbital periods of Jupiter.

N_dt_samples = 100
dt_samples = sim.particles[1].P * np.logspace(-3,0.,N_dt_samples)

Let’s run the simulations with the standard WH integrator first. We set safe_mode to 0 to speed up the calculation.

Emax_wh = np.zeros(N_dt_samples)
for i, dt in enumerate(dt_samples):
    sim_run = sim.copy() # make a copy of the simulation so we don't need to set a new one up every time
    sim_run.integrator = "whfast"
    sim_run.dt = dt
    sim.ri_whfast.safe_mode = False
    Emax_wh[i] = measure_energy(sim_run)
f,ax = plt.subplots(1,1,figsize=(7,5))
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("timestep [Jupiter years]")
ax.set_ylabel("relative energy error")
ax.plot(dt_samples/sim.particles[1].P,Emax_wh,label="WH")
ax.legend();
../_images/HighOrderSymplectic_10_0.png

We can see that the WH converges quadratically and reaches a precision of \(10^{-9}\) at a timestep of 0.001 Jupiter years. Let’s try the same with the WHCKL, SABACL4, and SABA(10,6,4) integrators.

Emax_saba4 = np.zeros(N_dt_samples)
for i, dt in enumerate(dt_samples):
    sim_run = sim.copy() # make a copy of the simulation so we don't need to set a new one up every time
    sim_run.integrator = "SABACL4"
    sim_run.ri_saba.safe_mode = False
    sim_run.dt = dt
    Emax_saba4[i] = measure_energy(sim_run)
Emax_saba1064 = np.zeros(N_dt_samples)
for i, dt in enumerate(dt_samples):
    sim_run = sim.copy()
    sim_run.integrator = "SABA(10,6,4)"
    sim_run.ri_whfast.safe_mode = False
    sim_run.dt = dt
    Emax_saba1064[i] = measure_energy(sim_run)
Emax_whckl = np.zeros(N_dt_samples)
for i, dt in enumerate(dt_samples):
    sim_run = sim.copy()
    sim_run.integrator = "WHCKL"
    sim_run.ri_whfast.safe_mode = False
    sim_run.dt = dt
    Emax_whckl[i] = measure_energy(sim_run)
/Users/rein/git/rebound/rebound/simulation.py:664: RuntimeWarning: WHFast convergence issue. Timestep is larger than at least one orbital period.
  warnings.warn(msg[1:], RuntimeWarning)

We get a warning message because the largest timestep we try is larger than the innermost orbital period. You would not want to use such a large timestep in an actual simulation, but we can ignore the message for this test.

f,ax = plt.subplots(1,1,figsize=(7,5))
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlabel("timestep [Jupiter years]")
ax.set_ylabel("relative energy error")
ax.plot(dt_samples/sim.particles[1].P,Emax_wh,label="WH")
ax.plot(dt_samples/sim.particles[1].P,Emax_saba4,label="SABACL4")
ax.plot(dt_samples/sim.particles[1].P,Emax_whckl,label="WHCKL")
ax.plot(dt_samples/sim.particles[1].P,Emax_saba1064,label="SABA(10,6,4)")
ax.legend();
../_images/HighOrderSymplectic_14_0.png

We can see that the higher order SABA4, WHCKL, and SABA(10,6,4) integrators are doing significantly better than the standard WH method. For very small timesteps the methods are limited by double floating point precision. The SABA methods appear to be slightly better than the WHCKL method above, but note that the SABA methods are also much slower than the WHCKL per timestep. Taking this into account, the WHCKL has a slight advantage for timesteps larger than 1% of the shortest orbital period. For extremely high accuracy, the SABA(10,6,4) method is faster.

Note that by default all integrators set safe_mode=False and keep_unsynchronized=False. Depending on your application, you might want to change these flags.