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?

[sharkdp]

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.

we don't have keyword arguments or default values, but we do have structs. there's not a lot of syntactic support for working with them so far, but maybe they could be a good option for implementing this options struct, and then we could provide functions that generate common sets of options for ease of use?

[Ryan-D-Gast]

we don't have keyword arguments or default values, but we do have structs. there's not a lot of syntactic support for working with them so far, but maybe they could be a good option for implementing this options struct, and then we could provide functions that generate common sets of options for ease of use?

Yes an option/setting struct would solve the problem. I did experiment with it but noticed like you said, that the syntactic support didn't seem easy enough for users to use (as in it felt like it made the API hard to work with). So if you add more syntactic support that would be amazing.

My preferred addition (if not already) is a way to initialize a struct of type Foo only listing items that you want to define as different then some predefined default. If that is already possible or added then I can make add a set high quality differential equation solvers for the standard library.

I can work on making do with what is available already and implement some solvers or wait tell some numbat language editions so I don't have to redo any work. Let me know what works for you.

[sharkdp]

My preferred addition (if not already) is a way to initialize a struct of type Foo only listing items that you want to define as different then some predefined default. If that is already possible or added then I can make add a set high quality differential equation solvers for the standard library.

That's not yet possible, but I definitely would like to have that, too.

I can work on making do with what is available already and implement some solvers or wait tell some numbat language editions so I don't have to redo any work. Let me know what works for you.

It's probably best to wait a bit. I am planning to merge this PR in it's current form soon.

[Ryan-D-Gast]

It's probably best to wait a bit. I am planning to merge this PR in it's current form soon.

I agree.

I think this language has a lot of potential for use in numerical computing. For now this implementation should be sufficient for most use cases. In the future when new language features are added, I will work on a more comprehensive implementation for the standard library (assuming that is ok with you).

I enjoyed our back and forth conversations on this matter. I look forward to working more with you in the future.