SciML / OrdinaryDiffEq.jl

High performance ordinary differential equation (ODE) and differential-algebraic equation (DAE) solvers, including neural ordinary differential equations (neural ODEs) and scientific machine learning (SciML)
https://diffeq.sciml.ai/latest/
Other
535 stars 202 forks source link

ImplicitEulerExtrapolation keywords #2371

Open ArnoStrouwen opened 1 month ago

ArnoStrouwen commented 1 month ago

I'm not sure this branch is doing anything: https://github.com/SciML/OrdinaryDiffEq.jl/blob/dc0e1e733fdd79a6c643cbd0c4f4be03c934cb91/lib/OrdinaryDiffEqExtrapolation/src/algorithms.jl#L50-L57 The lpad stuff seems overkill also.

ChrisRackauckas commented 1 month ago

Yeah that branch is trivial, seems like a mistake.

There's a few things with the extrapolation methods that are a bit odd right now. They are one of the best methods when tuned (for example https://docs.sciml.ai/SciMLBenchmarksOutput/stable/StiffODE/Pollution/) but the default tuning is pretty bad. @utkarsh530 can we take a second to work with @ArnoStrouwen here on documenting the tuning process, and change the default tuning?

We should include https://ieeexplore.ieee.org/abstract/document/9926357 in the documentation of the methods. The original resource is Hairer I and II but our implementation goes well beyond that.

These methods seem to only ever make sense if they are using the multithreading, as if you cannot use the multi threading they are less efficient than other methods for stiff and non-stiff. So when documenting them, documenting that threading = OrdinaryDiffEq.PolyesterThreads() on a multithreaded machine is a pretty crucial aspect.

Generally the go-to one should be ImplicitEulerBarycentricExtrapolation, I haven't seen the HW methods do good in any of the benchmarks for implicit methods. For explicit probably ExtrapolationMidpointHairerWanner.

Bumping the minimum order tends to improve multithreaded efficiency since it requires more calculations per step, and so it can take larger stepsizes and fill your CPU much better if you force it. init order should start higher than the minimum.

For implicit extrapolation, you only get better multithreading here if it's below the threshold to where BLAS LU factorization multithreads efficiently, which is around 200x200 matrices or so. So its sweet spot is somewhere around 20-500 ODEs. For explicit extrapolation, higher order Runge-Kutta methods like Verner methods just tend to do better. But both the implicit and explicit extrapolation methods are arbitrary order, and they also multithread better at higher orders, so there is a precision value at which they will simply be the best. For explicit extrapolation vs something like Vern9, that seems to be outside of Float128 tolerances. For implicit extrapolation, it's hard to tell how that matches up against a good Radau implementation, since we just got Radau9 and Hairer's radau doesn't have our improved linear algebra but it does really well with its adaptive order 5/9/13. We plan to complete the 5/9/13/beyond version this summer though, in which case it's really a toss up between what will be RadauIIA() vs ImplicitEulerBarycentricExtrapolation for the best method for extremely low tolerances in stiff equations.