gabrielgellner/DiffEQ.jl

Runge-Kutta method types

Closed this issue · 1 comments

Currently we closely follow the structure of the ODE.jl packages structure
which is extremely generic, allowing for inputs of any Butcher tableau with
a generic stepper, one for fixed step and one for adaptive. We have removed
the fixed steppers as it is of very little utility outside of very specific
use cases. In generally it increasingly feels like there is little practical
benefit for this generality. From the treatise by Hairer and Wanner (1993) it
is clear that Dopri5/Dop853 are a clear winners among the explicit Runge-Kutta
solvers, crushing RKFehlberg45 codes. The only exception might be the more
recent ode23 like codes from Shampine (which look like they might be more efficient for coarse tolerance and mild stiffness). Whatever the case we will move towards
only supporting these 2-3 solvers for the explicit Runge-Kutta family. To this
end some nice simplifications can be made:

  • Simplify the Butcher tableau type
    • The FSAL property should be a Bool field that is set at construction [Done]
    • The stages (S) of the tableau should be a field set at construction and
      not a type parameter [Done]

Once these changes are implemented the question comes on how to call the Dopri
codes. Currently we have specific types for the specific type in Dopri54.
No since we want to support both Dopri54 and Dopri853 the question becomes
whether it makes sense to use the current method or instead have a general
type Dopri that then has an argument for the order Dopri(func, order = 54).
And then for the regular case we just default to the best default solver
(Dopri853 I believe). The only issue with this is how to do dispatch in this
case as the solver method might need the order information for things like
the Hermite interpolation and how the step embedding works. The only solution
for this would be to either have manual dispatch with if statements. Or have
the order be type instances that can be delegated on, like
Dopri(func, order = dp54) where dp45 <: RKOrder or some such.

I really don't want to have any support for the tableaus anymore. It is a really cool piece of code, but it defaults the pragmatic view I have for this package. I am going to get rid of the tableau logic and just use basic arrays. The only trick will be to make sure I support having either Float32 or Float64.