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,
with initial conditions
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 thestep()function, this value represents the maximum timestep size in absolute value;callback: this is a callable with signaturebool (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
truethen 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 signaturevoid pre_hook(taylor_adaptive<double> &);
that will be invoked once before the first step is taken by the
propagate_for()andpropagate_until()functions.Step callbacks are passed by value to
propagate_for()andpropagate_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 thestep()function, this value represents the maximum timestep size in absolute value;callback: this is a callable with signaturebool (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
truethen 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 signaturevoid 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';
}