Customising the adaptive integrator

In the previous section we showed a few usage examples of the taylor_adaptive class using the default options. Here, we will show how the behaviour of the integrator can be customised in a variety of ways.

Error tolerance

As we mentioned earlier, by default the taylor_adaptive class uses an error tolerance equal to the machine epsilon of the floating-point type in use. E.g., when using double precision, the tolerance is set to \(\sim 2.2\times 10^{-16}\).

The tolerance value is used by the taylor_adaptive class to control the error arising from truncating the (infinite) Taylor series representing the solution of the ODE system. In other words, taylor_adaptive strives to ensure that the magnitude of the remainders of the Taylor series is not greater than the tolerance, either in an absolute or relative sense. Absolute error control mode is activated when all elements of the state vector have a magnitude less than 1, while relative error control mode is activated when at least one element of the state vector has a magnitude greater than 1.

In order to specify a non-default tolerance, the keyword argument tol can be used when constructing an integrator object:

import heyoka as hy

# Create the symbolic variables x and v.
x, v = hy.make_vars("x", "v")

# Create the integrator object.
ta = hy.taylor_adaptive(
                        # Definition of the ODE system:
                        # x' = v
                        # v' = -9.8 * sin(x)
                        sys = [(x, v),
                         (v, -9.8 * hy.sin(x))],
                        # Initial conditions for x and v.
                        state = [0.05, 0.025],
                        # Set the tolerance to 1e-9
                        tol = 1e-9)

ta

C++ datatype            : double
Tolerance               : 1e-09
High accuracy           : false
Compact mode            : false
Taylor order            : 12
Dimension               : 2
Time                    : 0
State                   : [0.05, 0.025]

The optimal Taylor order for a tolerance of \(10^{-9}\) is now 12 (instead of 20 for a tolerance of \(\sim 2.2\times 10^{-16}\)).

Integrating the system back and forth shows how the accuracy of the integration is reduced with respect to the default tolerance value:

# Integrate forth to t = 10 and then back to t = 0.
ta.propagate_until(t = 10.)
ta.propagate_until(t = 0.)

ta

C++ datatype            : double
Tolerance               : 1e-09
High accuracy           : false
Compact mode            : false
Taylor order            : 12
Dimension               : 2
Time                    : 0
State                   : [0.05000000000131285, 0.02499999999755865]

Compact mode

By default, the just-in-time compilation process of heyoka.py aims at maximising runtime performance over everything else. In practice, this means that heyoka.py generates a timestepper function in which there are no branches and where all loops have been fully unrolled.

This approach leads to highly optimised timestepper functions, but, on the other hand, it can result in long compilation times and high memory usage for large ODE systems. Thus, heyoka.py provides also a compact mode option in which code generation employs more traditional programming idioms that greatly reduce compilation time and memory usage. Compact mode results in a performance degradation of \(\lesssim 2\times\) with respect to the default code generation mode, but it renders heyoka.py usable with ODE systems consisting of thousands of terms.

Let’s try to quantify the performance difference in a concrete case. In this example, we first construct the ODE system corresponding to an N-body problem with 6 particles via the model.nbody() utility function:

# Create an nbody system with 6 particles.
sys = hy.model.nbody(n = 6)

Next, we create an initial state vector for our system. The contents of the vector do not matter at this stage:

# Create an initial state vector (6 values per body).
import numpy as np
sv = np.zeros(36)

Next, we time the creation of an integrator object in default code generation mode:

%time ta_default = hy.taylor_adaptive(sys, sv)

CPU times: user 3.2 s, sys: 36.5 ms, total: 3.24 s
Wall time: 3.24 s

Finally, we time the creation of the same integrator object in compact mode (which can be activated via the compact_mode keyword argument):

%time ta_default = hy.taylor_adaptive(sys, sv, compact_mode = True)

CPU times: user 211 ms, sys: 87 μs, total: 211 ms
Wall time: 207 ms

That is, in this specific example compact mode is more than 10 times faster than the default code generation mode when it comes to the construction of the integrator object. For larger ODE systems, the gap will be even wider.

High-accuracy mode

For long-term integrations at very low error tolerances, heyoka.py offers an opt-in high-accuracy mode. In high-accuracy mode, heyoka.py employs techniques that minimise the numerical errors arising from the use of finite-precision floating-point numbers, at the cost of a slight runtime performance degradation.

Currently, high-accuracy mode changes the way heyoka,py evaluates the Taylor polynomials used to update the state of the system at the end of an integration timestep. Specifically, polynomial evaluation via Horner’s rule is replaced by compensated summation, which prevents catastrophic cancellation issues and ultimately helps maintaining machine precision over very long integrations.

High-accuracy mode can be enabled via the high_accuracy keyword argument.