Notebook

Brouwer's law in the outer Solar System¶

In this tutorial we will show how heyoka.py can be used for long-term integrations of the outer Solar System accurate to machine precision.

Long-term integrations of the Solar System are often undertaken with symplectic integrators, which guarantee (from the point of view of the integration scheme) the conservation of dynamical invariants such as the total energy of the system. Because energy conservation is enforced by the integration scheme, in symplectic integrators the only source of error for the conservation of the total energy of the system derives from the use of approximate floating-point arithmetics. A result known as Brouwer's law establishes that the energy error deriving from the use of floating-point arithmetics cannot grow slower than $\sim \sqrt{t}$ (i.e., as a one-dimensional random walk). In other words, for any numerical integrator implemented on a real computer, the optimal behaviour (with respect to energy conservation in long-term integrations) is an error that grows with the square root of time.

Although heyoka.py is not a symplectic integrator, it is nevertheless able to achieve Brouwer's law, if properly configured. Specifically, in order to achieve Brouwer's law with heyoka.py we will need to:

lower the integration tolerance below machine precision, and
adopt techniques to reduce the numerical errors arising from the use of floating-point arithmetics.

The dynamical system¶

In this example we will study the dynamics the outer Solar System, that is, a 6-body problem consisting of the Sun, Jupiter, Saturn, Uranus, Neptune and Pluto, all considered as mutually-interacting point masses. We will adopt the Solar mass $M_\odot$ as unit of mass, the astronomical unit as unit of distance and the calendar year (365 days) as unit of time.

Let us begin with the definition of the numerical constants:

In [1]:

# Masses, from Sun to Pluto.
import numpy as np
masses = np.array([1.00000597682, 1 / 1047.355, 1 / 3501.6, 1 / 22869., 1 / 19314., 7.4074074e-09])

# The gravitational constant.
G = 0.01720209895 * 0.01720209895 * 365 * 365

Note how the Sun's mass is not exactly 1 because we included in it the mass of the terrestrial planets.

Next, we define a vector of cartesian initial conditions for the system. The numerical values are taken from this paper.

In [2]:

ic = [# Sun.
      -4.06428567034226e-3, -6.08813756435987e-3, -1.66162304225834e-6, +6.69048890636161e-6 * 365,
      -6.33922479583593e-6 * 365, -3.13202145590767e-9 * 365,
      # Jupiter.
      +3.40546614227466e+0, +3.62978190075864e+0, +3.42386261766577e-2, -5.59797969310664e-3 * 365,
      +5.51815399480116e-3 * 365, -2.66711392865591e-6 * 365,
      # Saturn.
      +6.60801554403466e+0, +6.38084674585064e+0, -1.36145963724542e-1, -4.17354020307064e-3 * 365,
      +3.99723751748116e-3 * 365, +1.67206320571441e-5 * 365,
      # Uranus.
      +1.11636331405597e+1, +1.60373479057256e+1, +3.61783279369958e-1, -3.25884806151064e-3 * 365,
      +2.06438412905916e-3 * 365, -2.17699042180559e-5 * 365,
      # Neptune.
      -3.01777243405203e+1, +1.91155314998064e+0, -1.53887595621042e-1, -2.17471785045538e-4 * 365,
      -3.11361111025884e-3 * 365, +3.58344705491441e-5 * 365,
      # Pluto.
      -2.13858977531573e+1, +3.20719104739886e+1, +2.49245689556096e+0, -1.76936577252484e-3 * 365,
      -2.06720938381724e-3 * 365, +6.58091931493844e-4 * 365]

We can now proceed to the definition of the dynamical equations. We will be using the make_nbody_sys() function, which sets up an ODE system corresponding to an N-body problem in cartesian coordinates:

In [3]:

import heyoka as hy
sys = hy.make_nbody_sys(6, masses = masses, Gconst = G)

The next step is the creation of the numerical integrator. We will be using a [batch integrator](<./Batch mode overview.ipynb>), which will allow us to substantially increase the floating-point throughput by integrating multiple sets of initial conditions at once. When creating the integrator, we will specify a tolerance of $10^{-18}$ (below machine precision) and we will activate high-accuracy mode. In high-accuracy mode, the integrator internally uses techniques (based on compensated summation and similar algorithms) to reduce the numerical errors arising from the use of floating-point arithmetics, at the price of a slight performance penalty.

In [4]:

# Multiplex the initial conditions to batches of 4 elements.
ic_batch = np.repeat(ic, 4).reshape(-1, 4)

# Create the integrator object, specifying a tolerance
# below machine precision and activating high-accuracy mode.
ta = hy.taylor_adaptive_batch(sys, ic_batch, high_accuracy = True, tol = 1e-18)

Integrating in parallel¶

In order to add statistical weight to our experiment, we will be integrating multiple sets of initial conditions at the same time. Each set of initial conditions will be slightly and randomly altered with respect to the numerical values introduced earlier, which will allow us to study the energy-conservation behaviour of the integrator using an ensemble of different but related problems.

As explained earlier, the use of a batch integrator already allows us to integrate multiple sets of initial conditions (4, in this case) at the same time. Additionally, we will concurrently run multiple batch integrators in different threads, in order to take full advantage of modern multi-core processors. The total integration time will be limited to 1 million years.

Let us take a look at the code:

In [5]:

# Define a logarithmic time grid over which
# the integrations will be performed.
t_grid = np.repeat(np.logspace(0, 6, 1000), 4).reshape(-1, 4)

# Multiplex the masses to batches of 4 elements.
masses_batch = np.repeat(masses, 4).reshape(-1, 4)

# A function for the computation of the total energy
# of the system from the state vector.
def energy(st):
        # Kinetic energy.
        vx = st[3::6]
        vy = st[4::6]
        vz = st[5::6]
        
        kin = np.sum(masses_batch * (vx**2 + vy**2 + vz**2) / 2, axis = 0)
        
        # Potential energy.
        pot = 0.
        for i in range(6):
            xi = st[i*6 + 0, :]
            yi = st[i*6 + 1, :]
            zi = st[i*6 + 2, :]
            
            for j in range(i+1, 6):
                xj = st[j*6 + 0, :]
                yj = st[j*6 + 1, :]
                zj = st[j*6 + 2, :]

                pot += -G * masses_batch[i] * masses_batch[j] / np.sqrt((xi - xj) * (xi - xj) + (yi - yj) * (yi - yj) + (zi - zj) * (zi - zj))

        return kin + pot

# The worker function that will be run in each thread.
def worker():
    from copy import deepcopy
    
    # Make a deep copy of the original integrator.
    ta_local = deepcopy(ta)
    
    # Randomly alter the initial conditions.
    new_state = ta_local.state + abs(ta_local.state) * np.random.uniform(-1e-12, 1e-12, ta_local.state.shape)
    
    # Determine the new centre of mass and its velocity.
    com_x = np.sum(new_state[0::6] * masses_batch, axis=0) / np.sum(masses_batch, axis=0)
    com_y = np.sum(new_state[1::6] * masses_batch, axis=0) / np.sum(masses_batch, axis=0)
    com_z = np.sum(new_state[2::6] * masses_batch, axis=0) / np.sum(masses_batch, axis=0)

    com_vx = np.sum(new_state[3::6] * masses_batch, axis=0) / np.sum(masses_batch, axis=0)
    com_vy = np.sum(new_state[4::6] * masses_batch, axis=0) / np.sum(masses_batch, axis=0)
    com_vz = np.sum(new_state[5::6] * masses_batch, axis=0) / np.sum(masses_batch, axis=0)
   
    # Recentre.
    new_state[0::6] -= com_x
    new_state[1::6] -= com_y
    new_state[2::6] -= com_z
    
    new_state[3::6] -= com_vx
    new_state[4::6] -= com_vy
    new_state[5::6] -= com_vz
   
    # Assign the new state.
    ta_local.state[:] = new_state
    
    # Compute the initial energy.
    E0 = energy(ta_local.state)

    # Integrate over the time grid.
    res = ta_local.propagate_grid(t_grid)
    
    # Check if any batch element produced an error.
    if any(oc[0] != hy.taylor_outcome.time_limit for oc in ta_local.propagate_res):
        raise RuntimeError("Integration failed: {}".format(ta_local.propagate_res))
    
    # Compute and return the relative energy error.
    return np.array([abs((E0 - energy(st)) / E0) for st in res])

The worker() function will be invoked concurrently from multiple threads of execution. It will first create a local copy of the integrator object, add some noise to the initial conditions, reset the centre of mass and then integrate the system for 1 million years.

Let us now run 16 batch integrations concurrently, for a total of $16\times 4 = 64$ sets of initial conditions.

NOTE: the following code will take a while to run (from a few seconds to a few minutes, depending on the CPU). In order to shorten the runtime, you can reduce the number of threads and/or the total integration time.

In [6]:

import concurrent.futures

# Run the integrations concurrently.
with concurrent.futures.ThreadPoolExecutor() as executor:
    futures = []
    for _ in range(16):
        futures.append(executor.submit(worker))

    # Gather the results.
    res = np.array([future.result() for future in concurrent.futures.as_completed(futures)])

Let us take a look at the shape of the result array:

In [7]:

res.shape

Out[7]:

(16, 1000, 4)

The first dimension refers to the 16 separate integrations we ran in parallel, the second dimension to the 1000 points in the time grid and the last dimension to the 4 elements in each batch. Let us re-arrange the array in order to facilitate further analysis:

In [8]:

res = res.transpose((1, 0, 2)).reshape((1000, -1))
res.shape

Out[8]:

(1000, 64)

Now the first dimension refers to the time grid points, and we have squashed into the second dimension the results of all integrations for each time point.

Results¶

We are now ready to plot the results of the integrations. For each time point, we will be plotting:

the relative energy error for all batch elements and parallel integrations,
the root mean square of the relative energy error across all batch elements and parallel integrations.

We will also add a dashed line representing Brouwer's law (i.e., $\sqrt{t}$) for comparison:

In [9]:

from matplotlib.pylab import plt
plt.rcParams["figure.figsize"] = (12,9)

# Use log scale on both axes.
plt.xscale('log')
plt.yscale('log')

# Plot all data points.
plt.plot(t_grid[:, 0], res, color='gray', marker='.', linestyle='None', alpha=.4, markersize=3.5)

# Plot the root mean square computed over all
# data points at each timestep.
plt.plot(t_grid[:, 0], np.sqrt(np.mean(res*res, axis = 1)), color='k', label="RMS")

# Plot sqrt(t).
plt.plot(t_grid[:, 0], 1.5e-16 * np.sqrt(t_grid[:, 0]), 'k--', label="$\sqrt{t}$ (Brouwer's law)")

plt.xlabel("Time (years)")
plt.ylabel("Rel. energy error")
plt.legend();

The plot clearly shows how the average energy error starts out around machine precision and then begins to grow following Brouwer's law. These results indicate that heyoka.py is able to optimally conserve the invariants of a dynamical system over long-term integrations.