Dense Output and t_eval

See examples of both these features in the “Demos” folder.

As of CyRK v0.10.0, support for dense output and accurate interpolation is implemented for all integrators when using pysolve_ivp and cysolve_ivp. Note, currently this feature is not available for nbsolve_ivp if that is a something you would like to see in the future please open a new issue on GitHub.

t_eval

Since CyRK uses an adaptive time stepping scheme it will only record results at time steps that have acceptable error (governed by the user-provided tolerances and the integration method). The resulting time domain (and associated y-values) will be non-uniform and may be more or less densely spaced than the user may wish. To correct this, a user may provide an array via the t_eval argument. CyRK will then compute its adaptive steps as usual, but will only save data for each t_eval point. If a t_eval point falls in-between the adaptive steps, then CyRK will use an interpolator (see next section) to determine a well-informed guess for the value.

Note

All elements of t_eval must fall between or at the border of th provided time_span.

Dense Output

In addition to the time domain, y-values, and any extra output (See Extra Output), CyRK can also build and store localized interpolator functions. These are constructed at each successful time step using the current state of the problem: derivatives, y-values, and time information. The interpolators can then be called with a desired time \(t\) and the resulting \(y\) will be produced, even if that specific \(t\) was never hit during integration (as long at \(t\) is within the provided time_span). These interpolators are what are used for estimating \(y\) values when t_eval is provided.

Note

The dense outputs are relatively large and must be heap allocated at each time step. They are also quite computationally to calculate. Leave dense_output=False unless required for your use case.

The performance hit when using only t_eval is much less because it only used one interpolator that is updated during the solution. No memory allocations or large objects are created. Note that the computation to rebuild the interpolator is still present so using t_eval will always be slower than not.

Interpolating Extra Outputs with Dense Output

As discussed in the Extra Output documentation, CyRK can capture additional outputs from the ODE solver process. These are non-dependent variables (non-\(y\) values) that are not used during integration error calculations but may be useful data for the user. This is triggered when num_extra is set to > 0 in either pysolve_ivp or cysolve_ivp. If dense_output is also set to True, then the final solution interpolators will also interpolate these extra outputs. This is done by making additional calls to the differential equation to determine what the values of the extra outputs are at each interpolated time step.

More details can be found in Extra Output documentation.

Using Dense Outputs

Example:

def cy_diffeq(dy, t, y):
    dy[0] = (1. - 0.01 * y[1]) * y[0]
    dy[1] = (0.02 * y[0] - 1.) * y[1]

    # If using extra output: set `num_extra=2`
    # dy[2] = (1. - 0.01 * y[1])
    # dy[3] = (0.02 * y[0] - 1.) 

import numpy as np
from CyRK import pysolve_ivp

initial_conds = np.asarray((20., 20.), dtype=np.complex128, order='C')
time_span = (0., 50.)
rtol = 1.0e-7
atol = 1.0e-8

result = \
    pysolve_ivp(cy_diffeq, time_span, initial_conds, method="RK45", dense_output=True,
                rtol=rtol, atol=atol,
                pass_dy_as_arg=True)
print("Was Integration was successful?", result.success)
print(result.message)
print("Size of solution: ", result.size)
import matplotlib.pyplot as plt
fig0, ax0 = plt.subplots()
ax0.plot(result.t, result.y[0], c='r')
ax0.plot(result.t, result.y[1], c='b')

# Call the result function to use the dense output interpolators.
new_time_array = np.linspace(5.0, 10.0, 250)
interpolated_y = result(new_time_array)

fig1, ax1 = plt.subplots()
ax1.plot(new_time_array, interpolated_y[0], c='r')
ax1.plot(new_time_array, interpolated_y[1], c='b')