# Why is RK45 used as the "default" method for non-stiff ODEs rather than a multistep one?

From what I read, the "default" go-to method for non-stiff ODEs is the Dormand-Prince Runge-Kutta pair; for instance, in Matlab docs, "Most of the time, ode45 should be the first solver you try".

However, from what I understand multistep methods should be able to achieve comparable accuracies (same order of convergence) with fewer function evaluations (one per step rather than seven per step with Dormand--Prince). I understand that changing step size is more difficult with multistep methods, but there are strategies to do that and I don't imagine it to be a major difficulty.

What is the logic behind Runge--Kutta methods being a better choice than multistep ones then? What compensates that 7-fold increase in the computational cost?

• There's also the issue that you need to start multistep schemes somehow. You don't have enough backward steps for high orders of accuracy in the first couple of time steps. Dec 10, 2021 at 18:24
• somewhat related in Space Exploration SE: Differences between numerical propagators (a few mentions of RK45)
– uhoh
Dec 11, 2021 at 1:34
• If correctly implemented, the function evaluations per step average to 6, as dopri5 is FSAL, first-same-as-last. Dec 11, 2021 at 7:43

First, let's establish that they are a good choice. The SciMLBenchmarks are probably the most comprehensive that there are as of right now for modern methods. This uses the vast number of methods available in DifferentialEquations.jl. Let's check the non-stiff ODEs:

Just to make sure, lets look at some benchmarks against MATLAB, Python, and R to see that these implementations are extremely optimized, outperforming those in MATLAB/Python/R by 50x+. So for fast implementations, across the board Runge-Kutta methods seem to do well on non-stiff ODEs and it's very rare when considering cost/error for multistep methods to really be competitive.

Why?

What compensates that 7-fold increase in the computational cost?

There are two major factors. One of them is stability. The explicit Adams-Bashforth methods are horrifically unstable, losing stability as they grow in order, and thus they are not very competitive in practice.

For this reason, optimized software for non-stiff ODEs which use multistep methods generally do not use the Adams-Bashforth methods, but instead they use the Adams-Bashforth-Moulton coefficients.

But wait, ABM is implicit, doesn't that mean it's for stiff equations? No, that's a very simplistic view. In fact, there's explicit integrators for stiff equations too, and so "implicit = stiff" is too simplistic of a view. Instead, implicit methods can be used for non-stiff equations by using less stable methods to solve the nonlinear system. I.e. you cheat a bit and don't use Newton (which is more stable) and instead use fixed-point iteration (with a predictor-corrector) form. However, once you've moved into this form, an update requires solving a fixed point iteration and thus "multistep methods take one f call to update" looks true from the equations but is not true in practice. Thus the calculus becomes a lot more fuzzy. Even then, the Adams-Bashforth-Moulton is still not very stable.

By 5th order you're at the stability like Euler's method... so while methods like CVODE with Adams coefficients, VCABM, ode113, etc. can technically all go to 13th order and take giant steps with few function iterations, in practice that doesn't happen often because of the stability bound.

For reference, look at the stability region for the Dormand-Prince 8/7 method:

This one isn't actually used in practice (contrary to common belief, the Dormand-Prince 8/7 method is not the method used in dop853, instead it's a different tableau that is only well-documented in the software itself, but I digress), but its stability region is still a lot larger. This means that if things start to become a little stability bound it can be taking larger steps for cost.

But what are all of those other numbers? That takes us to the second major factor: the leading coefficients of the truncation error. In contrast to what most mathematicians would hope, the purpose of numerical ODE solving is not to make $$\Delta t \rightarrow 0$$. In fact, it's the opposite: you want to make $$\Delta t$$ as large as possible (subject to hitting the error tolerance that you wished). Thus thinking about asymptotic $$\mathcal{O(\Delta t^n)}$$ is far too simplistic. People in the Runge-Kutta world noticed that if you optimize the leading coefficient of the truncation error, then the biggest error term is still pretty small even when $$\Delta t$$ is not zero, and this can allow for taking larger steps. Indeed, the Dormand-Prince method (the tableau behind ode45) is one of these "leading truncation error optimizing" methods. (In reality, it made an assumption and optimized under that constraint, while Tsitouras relaxed this constraint in 2011 and created a tableau that has 20% lower of a coefficient, which if you go back and look at the benchmarks you'll see that's Tsit5 which on average performs a little better than DP5, but that's details). But anyways, the multistep methods have pretty massive leading coefficient values compared to this which means that to hit error tolerance they can require a smaller $$\Delta t$$ and thus more work in the end.

Those two factors together are some of the big players. There's some other details that can matter too, but I think you can contribute most of the behavior to this.

Now, why are BDF methods like ode15s and CVODE okay for stiff equations? Won't the same leading truncation error argument apply? Yes, it does, but those methods are optimized on stiff equations because they take small steps. That's a whole different story for another time.

• I do not get why the type of the iterative solver matters. If the method is of a lesser order, just iterate longer, but you get the same solution, no? Or do you imply that only few iterations would be made by the "predictor-corrector"? But why can't be that done with Newton? Dec 11, 2021 at 7:21
• Thanks for this great insightful answer! One curiosity: what is the origin of the last image, the one with the red stability region of Dormand-Prince 8/7? Dec 11, 2021 at 8:12
• "But why can't be that done with Newton?" Good question. It can be done with Newton, but using a Newton method requires changing to a $\mathcal{O(n^3)}$ cost. In some cases, just one LU-factorization of the Jacobian in the Newton iterations can be more expensive than the entire ODE solve using other methods. The other methods are all linear in cost w.r.t. ODE system size. So Newton really is just for stiff ODEs for the stability unless the system is sufficiently small enough for this to not matter. Hence Adams codes using fixed point iteration by default. Dec 11, 2021 at 9:35
• Source for the Adams-Bashforth stability is this review book thing and the source for the DP87 stability is David Ketcheson's site Dec 11, 2021 at 9:37