Dense & continuous output#
Added in version 0.4.0.
One of the peculiar features of Taylor’s method is that it directly provides, via the Taylor series (4), dense (or continuous) output. That is, the Taylor series built by the integrator at each timestep can be used to compute the solution of the ODE system at any time within the timestep (and not only at the endpoint) via polynomial evaluation.
Because the construction of the Taylor series is part of the timestepping algorithm, support for dense output comes at essentially no extra cost. Additionally, because the dense output is computed via the Taylor series of the solution of the ODE system, its accuracy is guaranteed to respect the error tolerance set in the integrator.
Dense output can be used either from a low-level API, which gives direct access to the coefficients of the Taylor polynomials of the solution within a timestep, or from a higher-level API, which facilitates the common use case of using the coefficients of the Taylor polynomials to compute the continuous extension of the solution. If you are interested only in the latter, you can skip the next sections and jump directly to the continuous output section.
Dense output for the step() functions#
In order to illustrate how to use dense output in heyoka, we will keep things simple and go back to the simple pendulum:
// 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}};
Enabling dense output for the step() functions in heyoka is a two-step process.
The first step is to invoke one of the step() functions
passing an extra boolean parameter set to true:
// Integrate for a single timestep, and store
// the Taylor coefficients in the integrator.
ta.step(true);
The extra true function argument instructs the integrator
to record into an internal array the list of Taylor series
coefficients that were generated by the timestepping algorithm.
We can fetch a reference to the list of Taylor coefficients
via the get_tc() member function:
// Fetch the Taylor coefficients of order 0
// for x and v.
std::cout << "TC of order 0 for x: " << ta.get_tc()[0] << '\n';
std::cout << "TC of order 0 for v: " << ta.get_tc()[ta.get_order() + 1u] << '\n';
The Taylor coefficients are stored in row-major order,
each row referring to a different state variable. The number
of columns is the Taylor order of the integrator plus one.
Thus, index 0 in the array refers to the zero-order coefficient
for the x variable, while index ta.get_order() + 1 refers
to the zero-order coefficient for the v variable:
TC of order 0 for x: 0.05
TC of order 0 for v: 0.025
Indeed, the zero-order Taylor coefficients for the state variables are nothing but the initial conditions at the beginning of the timestep that was just taken.
Important
This last point is important and needs to be stressed again: the list of Taylor coefficients always refers to the last step taken and not to the next step that the integrator might take.
We are now ready to ask the integrator to compute the value of the solution at some arbitrary time. Let’s pick \(t = 0.1\), which is about halfway through the timestep that was just taken:
// Compute and print the dense output at t = 0.1.
auto &d_out = ta.update_d_output(0.1);
std::cout << "x(0.1) = " << d_out[0] << '\n';
std::cout << "y(0.1) = " << d_out[1] << '\n';
x(0.1) = 0.0500303
y(0.1) = -0.024398
The update_d_output() member function takes in input an absolute time coordinate
and returns a reference to an internal array that will contain the state of the system
at the specified time coordinate, as computed by the evaluation of the Taylor series.
update_d_output() can also be called with a time coordinate relative to the current
time by passing true as a second function argument.
Let’s now ask for the dense output at the very end of the timestep that was just taken, and let’s compare it to the current state vector:
// Compute the dense output at the end of the
// previous timestep and compare it to the current
// state vector.
ta.update_d_output(ta.get_time());
const auto &st = ta.get_state();
std::cout << "x rel. difference: " << (d_out[0] - st[0]) / st[0] << '\n';
std::cout << "v rel. difference: " << (d_out[1] - st[1]) / st[1] << '\n';
x rel. difference: 0
v rel. difference: -0
That is, as expected, the dense output at the end of the previous timestep matches the current state of the system to machine precision.
Before concluding, we need to highlight a couple of caveats regarding the use of dense output.
Important
First, it is the user’s responsibility to ensure that the array of Taylor
coefficients contains up-to-date values. In other words, the user needs to remember
to invoke the step() functions with the extra boolean argument set to true
before invoking update_d_output().
Failure to do so will result in update_d_output() producing incorrect values.
Important
Second, the accuracy of dense output is guaranteed to match the integrator’s
accuracy only if the time coordinate falls within the last step taken. Note that heyoka will
not prevent the invocation of update_d_output() with time coordinates outside the
guaranteed accuracy range - it is the user’s responsibility to be aware
that doing so will produce results whose accuracy does not match the integrator’s
error tolerance.
Dense output for the propagate_*() functions#
Added in version 0.8.0.
Dense output can be enabled also for the time-limited propagation functions
propagate_for() and propagate_until() via the boolean keyword argument
write_tc.
When write_tc is set to true, the propagate_*() functions
will internally invoke the step() function with the optional boolean
flag set to true, so that at the end of each timestep the Taylor coefficients
will be available. The Taylor coefficients can be used, e.g., inside the
step callback that can be passed to the propagate_*() functions.
Note that propagate_grid() always unconditionally writes the Taylor coefficients
at the end of each timestep, and thus using the write_tc argument is not necessary.
Continuous output#
Added in version 0.16.0.
The propagate_for() and propagate_until() functions can optionally return
a function object providing continuous output in the integration interval.
That is, this function object can be used to compute the solution at any time
within the time interval covered by propagate_for/until().
Let us see a concrete example of continuous_output in action:
// Reset time and state.
ta.get_state_data()[0] = 0.05;
ta.get_state_data()[1] = 0.025;
ta.set_time(0);
// Integrate up to t = 10, and request continuous output.
auto c_out = std::get<4>(ta.propagate_until(10., kw::c_output = true));
Continuous output can be requested via the c_output boolean keyword flag, and it is stored
as the fifth element of the tuple returned by propagate_for/until(). The continuous_output
function object
is wrapped in a std::optional (so that when continuous output is not requested, the
std::optional will be empty).
Let us try to print c_out to screen:
// Print to screen.
std::cout << *c_out << '\n';
Direction : forward
Time range: [0, 10)
N of steps: 48
The screen output informs us that the c_out object is capable of providing continuous output
in the \(\left[ 0, 10 \right)\) time interval, and that 48 steps were taken
during the integration.
Let us compute and print to screen the state of the system at a few time coordinates:
// Compute and print to screen the system state
// at different time coordinates.
for (auto tm : {0., 1.5, 4.3, 6.7, 8.9, 10.}) {
// Compute the state at tm.
(*c_out)(tm);
// Print it out:
std::cout << "time=" << tm << ", x=" << c_out->get_output()[0] //
<< ", v=" << c_out->get_output()[1] << '\n';
}
time=0, x=0.05, v=0.025
time=1.5, x=-0.0088572, v=0.156048
time=4.3, x=0.0375906, v=-0.106177
time=6.7, x=-0.0193535, v=-0.146456
time=8.9, x=-0.0424699, v=-0.0862923
time=10, x=0.0487397, v=0.0429423
The call operator of c_out takes a time value tm as input, and computes and writes to
an internal buffer the value of the state vector at tm. The get_output() member function
returns a reference to the internal buffer containing the state of the system at tm.
As we can see from the screen output, the state vector at \(t = 0\) corresponds to the initial conditions. The state vector at \(t=10\) (i.e., at the end of the integration interval) corresponds to the current state of the integrator object.
Continuous output is somewhat similar to propagate_grid(), in the sense that both allow to compute
the value of the solution at arbitrary time points. propagate_grid() is computationally more efficient,
but it requires to specify up-front the list of times at which the solution should be computed.
Continuous output, on the other hand, is not bound to a predetermined time grid, and it can thus be helpful
if time coordinates of interest can be identified only after the solution has been computed.
Before concluding, we need to highlight a couple of caveats regarding the
use of continuous_output.
Important
The continuous_output function object stores internally the time coordinate
and Taylor coefficients at the end of each step taken during the integration interval.
This means that the memory usage of a continuous_output object scales linearly
with the number of timesteps taken during the integration interval. Thus, for a sufficiently
long integration interval, the continuous_output object might end up exhausting
the available memory.
Important
Like for dense output, the accuracy of continuous_output is guaranteed to match the integrator’s
accuracy only if the time coordinate falls within the integration interval. Note that heyoka will
not prevent the use of a continuous_output object outside the
guaranteed accuracy range - it is the user’s responsibility to be aware
that doing so will produce results whose accuracy does not match the integrator’s
error tolerance.
Note that, like the other main classes in heyoka, continuous_output supports
serialisation.
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}};
// Integrate for a single timestep, and store
// the Taylor coefficients in the integrator.
ta.step(true);
// Fetch the Taylor coefficients of order 0
// for x and v.
std::cout << "TC of order 0 for x: " << ta.get_tc()[0] << '\n';
std::cout << "TC of order 0 for v: " << ta.get_tc()[ta.get_order() + 1u] << '\n';
// Compute and print the dense output at t = 0.1.
auto &d_out = ta.update_d_output(0.1);
std::cout << "x(0.1) = " << d_out[0] << '\n';
std::cout << "y(0.1) = " << d_out[1] << '\n';
// Compute the dense output at the end of the
// previous timestep and compare it to the current
// state vector.
ta.update_d_output(ta.get_time());
const auto &st = ta.get_state();
std::cout << "x rel. difference: " << (d_out[0] - st[0]) / st[0] << '\n';
std::cout << "v rel. difference: " << (d_out[1] - st[1]) / st[1] << '\n';
// Reset time and state.
ta.get_state_data()[0] = 0.05;
ta.get_state_data()[1] = 0.025;
ta.set_time(0);
// Integrate up to t = 10, and request continuous output.
auto c_out = std::get<4>(ta.propagate_until(10., kw::c_output = true));
// Print to screen.
std::cout << *c_out << '\n';
// Compute and print to screen the system state
// at different time coordinates.
for (auto tm : {0., 1.5, 4.3, 6.7, 8.9, 10.}) {
// Compute the state at tm.
(*c_out)(tm);
// Print it out:
std::cout << "time=" << tm << ", x=" << c_out->get_output()[0] //
<< ", v=" << c_out->get_output()[1] << '\n';
}
}