The adaptive integrator

The adaptive integrator#

The taylor_adaptive class provides an easy-to-use interface to heyoka’s main capabilities. Objects of this class can be constructed from a system of ODEs and a set of initial conditions (plus a number of optional configuration parameters with - hopefully - sensible defaults). Member functions are provided to propagate in time the state of the system, either step-by-step or by specifying time limits.

Let us see how we can use taylor_adaptive to integrate the ODE system of the simple pendulum,

\[\begin{split}\begin{cases} x^\prime = v \\ v^\prime = -9.8 \sin x \end{cases}\end{split}\]

with initial conditions

\[\begin{split}\begin{cases} x\left( 0 \right) = 0.05 \\ v\left( 0 \right) = 0.025 \end{cases}\end{split}\]

Construction#

Here’s the code:

    // Create the symbolic variables x and v.
    auto [x, v] = make_vars("x", "v");

    // Create the integrator object
    // in double precision.
    auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                      // x' = v
                                      // v' = -9.8 * sin(x)
                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                      // Initial conditions
                                      // for x and v.
                                      {0.05, 0.025}};

After creating the symbolic variables x and v, we construct an instance of taylor_adaptive<double> called ta. taylor_adaptive is a class template parametrised over the floating-point type which we want to use for the integration. In this case, we use the double type, meaning that the integration will be carried out in double precision.

As (mandatory) construction arguments, we pass in the system of differential equations using the syntax prime(x) = ..., and a set of initial conditions for x and v respectively.

Let us try to print to screen the integrator object:

    std::cout << ta << '\n';

This will produce the following output:

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

By default, the error tolerance of an adaptive integrator is set to the machine epsilon, which, for double, is \(\sim 2.2\times10^{-16}\). From this value, heyoka deduces an optimal Taylor order of 20, as indicated by the screen output. taylor_adaptive manages its own copy of the state vector and the time variable. Both the state vector and the time variable are updated automatically by the timestepping functions. Note also how, by default, the time variable is initially set to zero.

Single timestep#

Let us now try to perform a single integration timestep:

    // Perform a single step.
    auto [oc, h] = ta.step();

    // Print the outcome flag and the timestep used.
    std::cout << "Outcome : " << oc << '\n';
    std::cout << "Timestep: " << h << "\n\n";

    // Print again the integrator object to screen.
    std::cout << ta << '\n';

First, we invoke the step() member function, which returns a pair of values. The first value is a status flag indicating the outcome of the integration timestep, while the second value is the step size that was selected by heyoka in order to respect the desired error tolerance. We print both to screen, and we also print again the ta object in order to inspect how state and time have changed. The screen output will look something like this:

Outcome : taylor_outcome::success
Timestep: 0.216053

C++ datatype            : double
Tolerance               : 2.220446049250313e-16
High accuracy           : false
Compact mode            : false
Taylor order            : 20
Dimension               : 2
Time                    : 0.21605277478009474
State                   : [0.04399644836992638, -0.07844245547068798]

It is also possible to perform a single timestep backward in time via the step_backward() function:

    // Perform a step backward.
    std::tie(oc, h) = ta.step_backward();

    // Print the outcome flag and the timestep used.
    std::cout << "Outcome : " << oc << '\n';
    std::cout << "Timestep: " << h << "\n\n";

Outcome : taylor_outcome::success
Timestep: -0.213123

The step() function can also be called with an argument representing the maximum step size max_delta_t: if the adaptive timestep selected by heyoka is larger (in absolute value) than max_delta_t, then the timestep will be clamped to max_delta_t:

    // Perform a step forward in time, clamping
    // the timestep size to 0.01.
    std::tie(oc, h) = ta.step(0.01);

    // Print the outcome flag and the timestep used.
    std::cout << "Outcome : " << oc << '\n';
    std::cout << "Timestep: " << h << "\n\n";

    // Perform a step backward in time, clamping
    // the timestep size to 0.02.
    std::tie(oc, h) = ta.step(-0.02);

    // Print the outcome flag and the timestep used.
    std::cout << "Outcome : " << oc << '\n';
    std::cout << "Timestep: " << h << "\n\n";

Outcome : taylor_outcome::time_limit
Timestep: 0.01

Outcome : taylor_outcome::time_limit
Timestep: -0.02

Note that the integration outcome is now taylor_outcome::time_limit, instead of taylor_outcome::success.

Before moving on, we need to point out an important caveat when using the single step functions:

Warning

If the exact solution of the ODE system can be expressed as a polynomial function of time, the automatic timestep deduction algorithm may return a timestep of infinity. This is the case, for instance, when integrating the rectilinear motion of a free particle or the constant-gravity free-fall equation. In such cases, the step functions should be called with a finite max_delta_t argument, in order to clamp the timestep to a finite value.

Note that the propagate_*() functions (described below) are not affected by this issue.

Accessing state and time#

It is possible to read from and write to both the time variable and the state vector. The get_time()/set_time() functions can be used to access the time variable, while the get_state() and get_state_data() can be used to access the state data:

    // Print the current time.
    std::cout << "Current time: " << ta.get_time() << '\n';

    // Print out the current value of the x variable.
    std::cout << "Current x value: " << ta.get_state()[0] << "\n\n";

    // Reset the time and state to the initial values.
    ta.set_time(0.);
    ta.get_state_data()[0] = 0.05;
    ta.get_state_data()[1] = 0.025;

Note that get_state() returns a const reference to the std::vector holding the integrator state, while get_state_data() returns a naked pointer to the state data. Only get_state_data() can be used to mutate the state.

Time-limited propagation#

In addition to the step-by-step integration functions, taylor_adaptive also provides functions to propagate the state of the system for a specified amount of time. These functions are called propagate_for() and propagate_until(): the former integrates the system for a specified amount of time, the latter propagates the state up to a specified epoch.

Let us see a couple of usage examples:

    // Propagate for 5 time units.
    auto [status, min_h, max_h, nsteps, _1, _2] = ta.propagate_for(5.);

    std::cout << "Outcome      : " << status << '\n';
    std::cout << "Min. timestep: " << min_h << '\n';
    std::cout << "Max. timestep: " << max_h << '\n';
    std::cout << "Num. of steps: " << nsteps << '\n';
    std::cout << "Current time : " << ta.get_time() << "\n\n";

    // Propagate until t = 20.
    std::tie(status, min_h, max_h, nsteps, _1, _2) = ta.propagate_until(20.);

    std::cout << "Outcome      : " << status << '\n';
    std::cout << "Min. timestep: " << min_h << '\n';
    std::cout << "Max. timestep: " << max_h << '\n';
    std::cout << "Num. of steps: " << nsteps << '\n';
    std::cout << "Current time : " << ta.get_time() << "\n\n";

Outcome      : taylor_outcome::time_limit
Min. timestep: 0.202133
Max. timestep: 0.218136
Num. of steps: 24
Current time : 5

Outcome      : taylor_outcome::time_limit
Min. timestep: 0.202122
Max. timestep: 0.218139
Num. of steps: 72
Current time : 20

The time-limited propagation functions return a tuple of 6 values, which represent, respectively:

  • the outcome of the integration (which will usually be taylor_outcome::time_limit),

  • the minimum and maximum integration timesteps that were used in the propagation,

  • the total number of steps that were taken,

  • the continuous output function object, if requested (off by default),

  • a copy of the step callback (see below for an explanation). If no callback has been provided, an empty callback is returned.

The time-limited propagation functions can be used to propagate both forward and backward in time:

    // Propagate back to t = 0.
    std::tie(status, min_h, max_h, nsteps, _1, _2) = ta.propagate_until(0.);

    std::cout << "Outcome      : " << status << '\n';
    std::cout << "Min. timestep: " << min_h << '\n';
    std::cout << "Max. timestep: " << max_h << '\n';
    std::cout << "Num. of steps: " << nsteps << '\n';
    std::cout << "Current time : " << ta.get_time() << "\n\n";

    std::cout << ta << '\n';

Outcome      : taylor_outcome::time_limit
Min. timestep: 0.202078
Max. timestep: 0.21819
Num. of steps: 97
Current time : 0

C++ datatype            : double
Tolerance               : 2.220446049250313e-16
High accuracy           : false
Compact mode            : false
Taylor order            : 20
Dimension               : 2
Time                    : 0
State                   : [0.050000000000000044, 0.02499999999999999]

Note also that the time-limited propagation functions will stop integrating if a non-finite value is detected in the state vector at the end of the timestep. In such case, the outcome of the integration will be taylor_outcome::err_nf_state.

Added in version 0.7.0.

The propagate_for() and propagate_until() functions can be invoked with additional optional keyword arguments:

  • max_delta_t: similarly to the step() function, this value represents the maximum timestep size in absolute value;

  • callback: this is a callable with signature

    bool (taylor_adaptive<double> &);

    which will be invoked at the end of each timestep, with the integrator object as only argument. This is known as a step callback. If the callback returns true then the integration will continue after the invocation of the callback, otherwise the integration will be interrupted.

    Added in version 1.0.0.

    Optionally, a step callback can implement a pre_hook() member function with signature

    void pre_hook(taylor_adaptive<double> &);

    that will be invoked once before the first step is taken by the propagate_for() and propagate_until() functions.

    Step callbacks are passed by value to propagate_for() and propagate_until(), and they are returned as the sixth member of the return tuple;

  • c_output: a boolean flag that enables continuous output.

Propagation over a time grid#

Added in version 0.4.0.

Another way of propagating the state of a system in a taylor_adaptive integrator is over a time grid. In this mode, the integrator uses dense output to compute the state of the system over a grid of time coordinates provided by the user. If the grid is denser than the typical timestep size, this can be noticeably more efficient than repeatedly calling propagate_until() on the grid points, because propagating the system state via dense output is much faster than taking a full integration step.

Let us see a simple usage example:

    // Reset the time and state to the initial values.
    ta.set_time(0.);
    ta.get_state_data()[0] = 0.05;
    ta.get_state_data()[1] = 0.025;

    // Propagate over a time grid from 0 to 1
    // at regular intervals.
    auto out = ta.propagate_grid({0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0});

propagate_grid() takes in input a grid of time points represented as a std::vector, and returns a tuple of 6 values. The first 4 values are the same as in the other propagate_*() functions:

  • the outcome of the integration,

  • the minimum and maximum integration timesteps that were used in the propagation,

  • the total number of steps that were taken.

The fifth value returned by propagate_grid() is the step callback, if provided by the user (see below). Otherwise, an empty callback is returned.

The sixth value returned by propagate_grid() is a std::vector containing the state of the system at the time points in the grid. The state vectors are stored contiguously in row-major order:

    // Print the state at t = 0.4 (index 4 in the time grid).
    std::cout << "x(0.4) = " << std::get<5>(out)[2 * 4] << '\n';
    std::cout << "v(0.4) = " << std::get<5>(out)[2 * 4 + 1] << '\n';

x(0.4) = 0.0232578
v(0.4) = -0.14078

There are a few requirements on the time values in the grid:

  • they must all be finite,

  • they must be ordered monotonically,

  • the first value in the grid must be equal to the current time of the integrator.

Changed in version 4.0.0: The requirement on the first value of the time grid.

Added in version 0.7.0.

The propagate_grid() function can be invoked with additional optional keyword arguments:

  • max_delta_t: similarly to the step() function, this value represents the maximum timestep size in absolute value;

  • callback: this is a callable with signature

    bool (taylor_adaptive<double> &);

    which will be invoked at the end of each timestep, with the integrator object as only argument. This is known as a step callback. If the callback returns true then the integration will continue after the invocation of the callback, otherwise the integration will be interrupted.

    Added in version 1.0.0.

    Optionally, a step callback can implement a pre_hook() member function with signature

    void pre_hook(taylor_adaptive<double> &);

    that will be invoked once before the first step is taken by the propagate_grid() function.

    Step callbacks are passed by value to propagate_grid() and they are returned as the fifth member of the return tuple.

Full code listing#

#include <iostream>

#include <heyoka/heyoka.hpp>

using namespace heyoka;

int main()
{
    // Create the symbolic variables x and v.
    auto [x, v] = make_vars("x", "v");

    // Create the integrator object
    // in double precision.
    auto ta = taylor_adaptive<double>{// Definition of the ODE system:
                                      // x' = v
                                      // v' = -9.8 * sin(x)
                                      {prime(x) = v, prime(v) = -9.8 * sin(x)},
                                      // Initial conditions
                                      // for x and v.
                                      {0.05, 0.025}};

    // Let's print to screen the integrator object.
    std::cout << ta << '\n';

    // Perform a single step.
    auto [oc, h] = ta.step();

    // Print the outcome flag and the timestep used.
    std::cout << "Outcome : " << oc << '\n';
    std::cout << "Timestep: " << h << "\n\n";

    // Print again the integrator object to screen.
    std::cout << ta << '\n';

    // Perform a step backward.
    std::tie(oc, h) = ta.step_backward();

    // Print the outcome flag and the timestep used.
    std::cout << "Outcome : " << oc << '\n';
    std::cout << "Timestep: " << h << "\n\n";

    // Perform a step forward in time, clamping
    // the timestep size to 0.01.
    std::tie(oc, h) = ta.step(0.01);

    // Print the outcome flag and the timestep used.
    std::cout << "Outcome : " << oc << '\n';
    std::cout << "Timestep: " << h << "\n\n";

    // Perform a step backward in time, clamping
    // the timestep size to 0.02.
    std::tie(oc, h) = ta.step(-0.02);

    // Print the outcome flag and the timestep used.
    std::cout << "Outcome : " << oc << '\n';
    std::cout << "Timestep: " << h << "\n\n";

    // Print the current time.
    std::cout << "Current time: " << ta.get_time() << '\n';

    // Print out the current value of the x variable.
    std::cout << "Current x value: " << ta.get_state()[0] << "\n\n";

    // Reset the time and state to the initial values.
    ta.set_time(0.);
    ta.get_state_data()[0] = 0.05;
    ta.get_state_data()[1] = 0.025;

    // Propagate for 5 time units.
    auto [status, min_h, max_h, nsteps, _1, _2] = ta.propagate_for(5.);

    std::cout << "Outcome      : " << status << '\n';
    std::cout << "Min. timestep: " << min_h << '\n';
    std::cout << "Max. timestep: " << max_h << '\n';
    std::cout << "Num. of steps: " << nsteps << '\n';
    std::cout << "Current time : " << ta.get_time() << "\n\n";

    // Propagate until t = 20.
    std::tie(status, min_h, max_h, nsteps, _1, _2) = ta.propagate_until(20.);

    std::cout << "Outcome      : " << status << '\n';
    std::cout << "Min. timestep: " << min_h << '\n';
    std::cout << "Max. timestep: " << max_h << '\n';
    std::cout << "Num. of steps: " << nsteps << '\n';
    std::cout << "Current time : " << ta.get_time() << "\n\n";

    // Propagate back to t = 0.
    std::tie(status, min_h, max_h, nsteps, _1, _2) = ta.propagate_until(0.);

    std::cout << "Outcome      : " << status << '\n';
    std::cout << "Min. timestep: " << min_h << '\n';
    std::cout << "Max. timestep: " << max_h << '\n';
    std::cout << "Num. of steps: " << nsteps << '\n';
    std::cout << "Current time : " << ta.get_time() << "\n\n";

    std::cout << ta << '\n';

    // Reset the time and state to the initial values.
    ta.set_time(0.);
    ta.get_state_data()[0] = 0.05;
    ta.get_state_data()[1] = 0.025;

    // Propagate over a time grid from 0 to 1
    // at regular intervals.
    auto out = ta.propagate_grid({0., 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0});

    // Print the state at t = 0.4 (index 4 in the time grid).
    std::cout << "x(0.4) = " << std::get<5>(out)[2 * 4] << '\n';
    std::cout << "v(0.4) = " << std::get<5>(out)[2 * 4 + 1] << '\n';
}