3 ms·
I recently wrote an ODE IVP integrator in C#. In the process I read the code in Numerical Recipes, Scipy and DifferentialEquations.jl. That combined with my o
by lamontcg 3y ago
I recently wrote an ODE IVP integrator in C#. In the process I read the code in Numerical Recipes, Scipy and DifferentialEquations.jl. That combined with my own house style and some gyrations for allocation-free code resulted in the final product. A lot of what I took was just how to abstract the code into a few different layers to make it a bit more pluggable so that now I've got swappable BS5 and DP5 integrators. That layer structure came mostly from scipy. Some of the details of DP5 and the PI stepsize control came from NR. Not a whole lot came from the julia code since it is honestly kind of a mess due to what needs to happen for performance there. I did dig into it enough to find out that they're just using bisection with left/right endpoint control for events, so I used that which was much simpler than trying to sort out what endpoint control with Brent's method looks like. Oh I took some tips of the overall shape of the events API from the julia code so that the events get the integrator itself injected into it, so events can tweak the integrator shape. Still got some details left that I need to clean up and flesh out the events API, but it is way better code than the slightly OOP'd wikipedia-quality code that I was using. It should be flexible enough to drop other RK methods into it and even drop non-RK IVP solvers into it.
Reasonably happy with the result and it is a lot better for having looked at other code and not continued to pursue NIH syndrome.