Add adaptive Runge-Kutta ODE solver (Dormand-Prince)

[Ryan-D-Gast]

Summary

Adds dsolve_adaptive_runge_kutta, an adaptive ODE solver using the Dormand-Prince 5(4) method.

Changes

• New dsolve_adaptive_runge_kutta(f, x_0, x_e, y_0, atol, rtol) function
• Hairer & Wanner initial step size estimator
• Updated documentation and tests

Why add this?

The existing dsolve_runge_kutta requires users to guess an appropriate step count. Too few steps gives poor accuracy; too many wastes computation. This is especially problematic for problems where the solution changes rapidly in some regions and slowly in others.

An adaptive solver automatically adjusts step sizes to maintain accuracy while minimizing work. Users specify error tolerances (atol, rtol) rather than step counts, making it both easier to use and more reliable.

Dormand-Prince is the standard choice for adaptive RK methods (used by MATLAB's ode45, SciPy's solve_ivp default, and Julia's OrdinaryDiffEq.jl).

Side Notes

1. Numbat currently lacks optional/default parameters and keyword arguments, making it difficult to expose solver configuration (safety factor, min/max step scaling, etc.) in a clean API. The current implementation uses SciPy's defaults. If these language features are added in the future, the solver API could be extended to allow user customization.
2. Second if the language supports it which I think it does in the future I could update this function to include a "solout" function which would allow for user configurable control of the output between steps (useful for getting even points and/or recording key events).

[sharkdp]

Thank you very much for your contribution!

I'm certainly not an expert on numerical algorithms. Would it make sense to replace my initial attempt with your solution?

[Ryan-D-Gast]

Hey,

If you look at my github profile you can see I've done a lot of work with numerical methods so maybe not an expert but I have a lot of experience.

The Dormand Prince algorithm is significantly more efficient and accurate and better for most use cases (hence why its default for SciPy). But I still think the rk4 algorithm is good to have as its more efficient for problems that do not require a high level of precision such as circuits and guidance navigate and control.

So the short answer no I think having both is ideal.

[Ryan-D-Gast]

In an ideal world I would create a function called solve_ivp aka just having a very similar API to SciPy's solve_ivp which would look like:

solve_ivp(fun, x0, xf, y0, method='RK45', solout=(some default output function which just logs the steps), **options)

Where **options is a struct or in the case of python a key word argument placeholder for changing settings of the solver such as the relative or absolute tolerance.

I would be happy to implement this which I've done before in rust. The problem is I think the languages features currently makes this near impossible.

Would it be possible or a good idea to implement it in Rust then expose it as a numbat function?