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';
    }
}

Dense & continuous output

Contents

Dense & continuous output#

Dense output for the step() functions#

Dense output for the propagate_*() functions#

Continuous output#

Full code listing#

Dense output for the `step()` functions#

Dense output for the `propagate_*()` functions#