4
$\begingroup$

How well do explicit Runge-Kutta "tableau" methods compare to the state of the art ODE solvers and when do they fail? I've been reading Butcher's ODE book and he does a good job at introducing tableaus to represent different Runge-Kutta methods. This includes methods with built-in error estimates for adaptive steps sizes like Fehlberg, Vergner, or Dorman-Prince.

What I don't understand is how these methods compare to the current state-of-the-art in ODE solvers. Are these methods ever really used or are they simply discussed for historical reasons? If they are used, what's their general limitation or domain of use?

As a minor aside, I noticed that DifferentialEquations.jl appeared to implement just about every known tableau method. That said, they don't appear to be the default methods used, which also caused me to question whether or not they're still of modern use.

Improve this question
$\endgroup$
6
  • $\begingroup$ My experience with higher-order RK methods for solid mechanics has been that the implementation complexity and extra computational cost can rarely be justified except for special cases. That may not be the case when you're solving, e.g., chemical rate equations. $\endgroup$
    – Biswajit Banerjee
    Sep 27, 2018 at 22:17
  • $\begingroup$ Can you explain what you mean by "state-of-the-art"? Runge-Kutta methods are one of the main classes of time integrators and many "state-of-the-art" methods are Runge-Kutta methods, so your question doesn't make sense. $\endgroup$
    – David Ketcheson
    Sep 28, 2018 at 4:40
  • 1
    $\begingroup$ @David Ketcheson For someone unfamiliar with ODEs, it's not a priori clear what classes retain their power in a modern setting. Wouldn't it be possible that variable-order Adams methods outperform RK in all situations, and that RK methods are all obsolete and just presented for historical reasons? (Benchmarks show this is false.) For example, in the comments here we see someone moving from one bad option to another bad option with a total misunderstanding of what to look for. $\endgroup$
    – procrastinating on actual work
    Sep 28, 2018 at 4:55
  • $\begingroup$ @procrastinatingonactualwork That's fine. The problem with your question, though, is that you draw a false dichotomy between Runge-Kutta methods and state-of-the-art methods. Also, the Butcher tableau is just a way of representing a Runge-Kutta method. Calling them Runge-Kutta tableau methods is confusing since somebody might understand that you don't simply mean Runge-Kutta methods. $\endgroup$
    – David Ketcheson
    Sep 28, 2018 at 6:59
  • $\begingroup$ @DavidKetcheson Sorry for the confusion. I'm unfamiliar with these methods and it's not my intention to minimize their impact. Butcher made a nice figure on page 90 of his book of different strategies for ODE solvers such as using past values (multistep), more calculations per step (Runge-Kutta), y derivatives (Taylor series), or f derivatives (Rosenbrock). He also introduced Almost Runge-Kutta (ARK) that appeared to be a hybrid method. I'm trying to understand whether taking RK coefficients and implementing them directly is 'good enough' or if there's some trick or hybrid that's better. $\endgroup$
    – wyer33
    Sep 28, 2018 at 13:33

3 Answers 3

Reset to default
8
$\begingroup$

I think the confusion here is what exactly is a "tableau method" since it's used in our docs in a very specific way. A "tableau method" in DifferentialEquations.jl parlance is a method which is implemented by explicitly building the arrays for the tableaus and performing loops on said arrays. The tableaus are stored in DiffEqDevTools.jl in the highest precision published and transform at compile-time to the right precision before solving. Here is Tsitouras' 5th order method:

https://github.com/JuliaDiffEq/DiffEqDevTools.jl/blob/master/src/ode_tableaus.jl#L866-L920

You can directly use these with the ExplicitRK method. Here is its internal perform_step function:

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/perform_step/explicit_rk_perform_step.jl#L72-L119

However, as another post mentioned, when we first did this we did a literature review of all of the tableaus we could find, implemented them, tested them, looked at stability, etc. (that 8000 line page of tableaus :) ) and then specialized some of the algorithms. For example, for Tsitoras' method we pulled it apart into a stack-allocated tableau

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/tableaus/low_order_rk_tableaus.jl#L497-L601

specialized the equations on the form of the constants

https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/perform_step/low_order_rk_perform_step.jl#L609-L656

added a special interpolant to it, modified some default magic numbers for adaptivity, etc. to specifically optimize this algorithm. We also have a trait system which for example hard codes specific constants about the stability region necessary for performing stiffness detection and automatic switching to implicit methods. And we did that for a few other RK methods. A lot of these optimizations only truly matter if f is inexpensive, but since that shows up in a lot of cases (parameter estimation and searches, analysis of chaotic systems) we took the time to do this, tested it, and performance-wise would recommend these implementations over the pure tableau-based form. Note that it is mathematically the same method, but this formulation does not require loops over constants and instead is in a form that LLVM optimizes very well. When the equation is small, Julia will in fact in line the user's f derivative function and the solver will construct something entirely different and efficient.

But that's just optimizations and just a DifferentialEquations.jl thing. Explicit RK methods, implicit RK methods, Rosenbrock methods, etc. can all be discussed as methods in a tableau form, and these are all quite good methods. The Tsitorious method which was shown above is actually one of the core defaults of the ecosystem, and a similar method is the Dormand-Prince method which you might know of as ode45. These are still quite state-of-the-art (at least for non-stiff equations, and when paired with PI-adaptivity and their special interpolants). These are the methods you see recommended everywhere, along with a few non-tableau methods like multistep Adams and BDF. Due to the benchmarks we have we recommend the hand-tuned versions, but it's stepping behavior are mathematically the same as if you looped through arrays of constants (sans the specializations on constants and interpolations).

For your aside, our test showed that not every tableau is worth it. Some just aren't efficient. Over time, people have gotten better at choosing coefficients which reduce the leading truncation error term while keeping the stability region the same size. The older methods then take similar steps but result in a solution with more error. We keep them there for research purposes, and also it's a nice archive of the literature since some of those are hard to find.

Improve this answer
$\endgroup$
0
4
$\begingroup$

Runge-Kutta methods are arguably the most widely used methods for initial value ODEs (including also PDE semi-discretizations). There are particular types of Runge-Kutta methods that are applicable to virtually every class of initial value problem. They're certainly not just a historical artifact.

Improve this answer
$\endgroup$
3
$\begingroup$

DifferentialEquations.jl has the best survey of current solvers. It goes through every option and benchmarks them. Modern RK tableaus are the state-of-the-art in most non-stiff situations. Multistep methods are good with expensive non-stiff functions. For stiff problems, it varies. See:

http://diffeq.sciml.ai/latest/solvers/ode_solve.html

http://devdocs.juliadiffeq.org/latest/internals/notes_on_algorithms.html

The two documents occasionally contradict each other, such as recommending different algorithms at order 7 and 8, so you'll have to ask Chris Rackauckas about that (see Chris's comment below for the reason). On the second link, there's an implicit understanding that higher order methods are useful when the desired tolerance shrinks, but also come with speed drawbacks at looser tolerances.

You'll also want to learn about stepsize control since many of the best algorithms rely on it and leave it as a separate problem. DiffEq has this short description but I think the original paper is more understandable.

Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ The "notes on algorithms" was a literature review which was done without the modern Verner algorithms, so it's all but the modern Verner. Modern Verner perform the best for orders 6-9. Found here people.math.sfu.ca/~jverner $\endgroup$
    – Chris Rackauckas
    Sep 28, 2018 at 16:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.