t_eval of CyRK behaves different to solve_ivp

[fvollmer]

I just figured out that t_eval isn't behaving in the same way as scipy.integrate.solve_ivp:

solve_ivp solves the ode for the points in t_eval

Times at which to store the computed solution, must be sorted and lie within t_span. If None (default), use points selected by the solver.

CyRK interpolates the solution t_eval

Both solvers uses an adaptive time stepping protocol based on the recent error at each step. This results in a final non-uniform time domain of variable size. If the user would like the results at specific time steps then they can provide a np.ndarray array at the desired steps via t_eval. The solver will then interpolate the results to fit this array.

In my kind of special case CyRK is actually slower than solve_ivp.

[jrenaud90]

Hi, @fvollmer! Yes this is a missing feature of CyRK. It would not be too difficult to implement: solve_ivp uses a dense output which estimates the value of your problem inbetween solution steps (sorta like a faster, more accurate interpolation).

See here for the main implementation in SciPy and here for the RK specific implementation.

• Update RK constants to include the versions used for dense output.
• For nbrk: Create a python class that largely mimics the dense output class used in SciPy. Add this as an optional return to nbrk_ode.
• For cyrk and CySolver: Create a cython class that does the same. For cyrk_ode this can be an additional output. For CySolver this class could be partially instantiated early on and then stored in the CySolver class as <instance_name>.dense and used like <instance_name>.dense(new_x) -> new_y.

I will mark this as a todo for a future update (PRs are welcome if you or others want to work on this!)

In the mean time feel free to share a simple working example of your code. I would be curious why the interpolation is so severely impacting performance. If you are not using the CySolver approach then that would be a way to greatly improve performance.

Cheers!