Adaptive RK4 Orbits

We’ll use an implementation of the orbit problem with RK4 written as a module. The source (in both C++ and python) is here:

zingale/computational_astrophysics

import numpy as np
import matplotlib.pyplot as plt

import orbit_adaptive as orbit

Let’s first look at the source

#%cat orbit_adaptive.py

Highly eccentric orbit

Let’s create an orbit with a tolerance of \(10^{-5}\), and the a very high eccentricity (\(a = 1\), \(e = 0.95\))

o = orbit.OrbitsRK4(1.0, 0.95)
o.integrate(0.05, 1.e-5, o.kepler_period())

How many points did it take?

o.npts()

92

Let’s plot it

fig = o.plot(points=True)
fig.set_size_inches(5, 12)

Clearly we see that the size of the timestep is much larger at aphelion as compared to perihelion.

Let’s create a version with a fixed timestep, with a timestep of \(0.001~\mathrm{yr}\). We can do this by passing in a negative error.

o_fixed = orbit.OrbitsRK4(1.0, 0.95)
o_fixed.integrate(0.001, -1, o_fixed.kepler_period())
fig = o_fixed.plot()
fig.set_size_inches(5, 12)

print(o_fixed.npts())

1001

We need a really small fixed step size to integrate this at all reasonably

Timestep evolution

Let’s look at the adaptive version some more. First, let’s look at the timestep evolution.

ts = np.array(o.time)

dt = ts[1:] - ts[:-1]
tc = 0.5 * (ts[1:] + ts[:-1])

fig, ax = plt.subplots()
ax.plot(tc, dt)
ax.set_xlabel("t [yr]")
ax.set_ylabel("dt [yr]")
ax.set_yscale("log")

Notice that the timestep changes by ~ 3 orders of magnitude over the evolution.

Energy conservation

What about energy conservation?

e = []
ts = o.time
for n in range(o.npts()):
    e.append(o.energy(n))
    
fig, ax = plt.subplots()
ax.plot(ts, e/e[0])
ax.set_xlabel("time [yr]")
ax.set_ylabel("E/E(t=0)")
ax.ticklabel_format(useOffset=False)

Now let’s do 10 orbits

o = orbit.OrbitsRK4(1.0, 0.95)
o.integrate(0.05, 1.e-5, 20*o.kepler_period())

e = []
ts = o.time
for n in range(o.npts()):
    e.append(o.energy(n))
    
fig, ax = plt.subplots()
ax.plot(ts, e/e[0])
ax.set_xlabel("time [yr]")
ax.set_ylabel("E/E(t=0)")
ax.ticklabel_format(useOffset=False)

Caution

We can see that RK4 does not conserve energy!

There is a steady decrease in the total energy over 20 orbits. If we wanted to evolve for millions of years, this would certainly be a problem. Making the tolerance tighter certainly will help, but it will also make things a lot more expensive.

The problem is that 4th order Runge-Kutta does not know anything about energy conservation.