So there is a ton to say about this, and we will actually be putting a paper out that tries to summarize it a bit, but let me narrow it down to something that can be put into a quick StackOverflow post. I will make one statement really early and keep repeating it: you cannot untangle the efficiency of a method from the efficiency of a software. The details of software implementations are what actually matter in this field.
Essentially, the state of the art methods completely depend on the problem that you're looking at. DiffEqBenchmarks.jl has quite a few benchmarks of some of the latest software (includes DifferentialEquations.jl, SUNDIALS, Hairer's FORTRAN stuff, Shampine's FORTRAN stuff, LSODA, etc.). It's based on Hairer's benchmarks from Hairer II but the DiffEqBenchmarks are more up-to-date with the software, and explores many of the methods mentioned here more thoroughly (while still including Hairer's FORTRAN methods). The full list of algorithms that are available to benchmark through this can be found in the DifferentialEquations.jl documentation, which covers most known algorithms at this point (and getting more complete after this summer). Not all of the methods in the documentation are showing up in every plot, but you can download the benchmarks and add/subtract things plus change around options.
Although you want to decouple "software" from the "methods", the actual software implementation of many details matters as much if not more than the actual method that is used. So you can see in many of the benchmarks on stiff problems that the three different BDF implementations (ddebdf in Fortran, LSODA, CVODE_BDF) have very different efficiencies on different problems given how they choose time steps, re-use Jacobians, etc. Those more "engineering" aspects (especially CVODE's Jacobian re-use structure) is what makes a software efficient or inefficient. And as you can see from the benchmarks, changing around options matters a ton as well. For example, SUNDIALS ARKODE fails to solve most of the stiff benchmarks we threw at it without changing some options (This is noted in the SUNDIALS examples as well, so it's not just an error on our end)! So keep that in mind while talking about methods and whenever reading a paper that talks about the efficiency.
Given all of those caveats, high order Rosenbrock methods seem to dominate when you have sufficiently small systems because they can take large time steps and don't require the stability of Newton iterations (which places the real stability bound on implicit methods). While previous literature mentioned that Jacobian accuracy issues would prevent Rosenbrock methods from converging well at high order, modern automatic differentiation techniques like those in ForwardDiff.jl circumvent a lot of the traditional issues with Jacobian accuracy of numerical differentiation which has breathed new life into this field (and is one of the reasons why Julia is being used in this domain). Parallelization of BLAS seems to only get in the way at this size (Jacobian <50x50, see, engineering details matter!), so I suspect there is some gains to be had looking at parallel Rosenbrock methods (which allow for multiple lu-factorizations to be done in parallel by creating a DIRK-like Rosenbrock), there aren't any real software implementations to test this idea.We will have software coming out soon enough to start benchmarking with, but the literature has only created methods for dual core parallel Rosenbrock methods (that I know of), so methods specialized for modern 4-8 core CPUs of order 4 or 5 would be an interesting topic. Parallelized impicit extrapolation methods could also be an interesting method in this domain.
One caveat I'll throw in here is that, as you ask for lower and lower tolerances (say around 1e-8), you'll always see Hairer's Radau implementation become the most efficient. It doesn't scale well to larger ODE sizes either, since it uses a Jacobian that is massively larger than the other methods as a FIRK method, but its higher order does become a factor that leads to impressive efficiency when high accuracy is required.
As the system size increasing to >100 ODEs, other factors start asymptotically mattering more. Essentially, the cost of computing a Jacobian and factorizing the W matrix $W = I - \gamma f'$ for some $\gamma \in \mathbb{R}$ become the dominating factors in any semi-implicit or implicit method (so any implicit RK, implicit multistep, Rosenbrock). As this occurs, the ability to re-use a Jacobian that you've already computed, or re-use a factorized Jacobian you've already created, becomes very important. Thus, since the accuracy of a standard Rosenbrock method depends on the accuracy of its Jacobian, Rosenbrock methods begin to lag behind because they need to compute and invert a new Jacobian every step. Implicit methods do not have this issue since the Jacobian is only used as a line search for solving the implicit equations, and thus does not actually effect accuracy. "Bad Jacobians" do cause more linear solves to have to be done in order for the Newton method to converge, but a linear solve with a pre-factorized matrix is asymptotically trivial in comparison to the Jacobian calculation or the inversion. Backsubstitution is $O(n^2)$, Jacobian calculations are $O(n^2)$, while inversion is $O(n^3)$. That fact makes people generally consider only the latter to be an issue, actual timing and profiling proves that is an oversimplification since calculating the Jacobian requires performing the calculations in $f$, the ODE derivative function, while the inversion is pure linear algebra (matrix coloring and analytical forms of sparse Jacobians reduce the complexity, but still require the extra calculations of $f$ in some form). Thus at this stage you have to segment the problems a little bit more based on how costly it is to compute $f$, since even if $N=100$ you could be given a really really slow $f$ calculation so $f$ calls may dominate time over the matrix inversions (this is generally the case in MATLAB or Python ODE programs).
So let's assume that you have an $f$-dominated semi-large ODE system, where we define semi-large to be a size that is large enough to care about these more asymptotic factors but small enough that dense/sparse factorizations still fit into memory. What we have found works well in these cases are things like SDIRK methods since they can stably step (L-stability and B-stability) with much higher order than BDF while re-using Jacobians to some extent. What happens is this minimizes $f$ calculations (though currently does require a few more inversions, something we are investigating, again an engineering "post-math" challenge) because BDF becomes efficient by doing a lot of Jacobian re-use and slamming really small steps (ha, bet you didn't see that coming) with the same inverted $W$ for multiple steps in a row. We need to formalize this into DiffEqBenchmarks.jl some more, but one good source on this is this plot from a PR which is looking at a stiff 1000 ODE model with costly f evals. In this case, KenCarp4 seems to be one of the better methods, though we are investigating alternative (E)SDIRK tableaus but haven't found a better one yet (none of the 5th order seem to do as well, a fact noted by Kennedy and Carpenter in their paper on ESDIRK methods as well). An alternative strategy which may prove useful here are RosW methods, which are Rosenbrock methods that do not lose order when the Jacobian is inaccurate. Comprehensive testing of these methods in this domain is set to occur during this summer. A lot of the efficiency will be dependent on developing efficient Jacobian re-use strategies though.
When finally getting to the domain that is actually dominated by matrix inversions, BDF seems to do well for weirder reasons than one might expect. Not BDF, but specifically the VODE and CVODE line of implementations. The reason is because they take very small steps in comparison to other stiff solvers, in a way that constrains their dt changes, which helps keep the $W$ matrix the same ($\gamma$ is always proportional to dt), and decreases the number of inversions which are required. This coupled with their special Fixed Leading Coefficient (FLC) form lets them re-use already inverted matrices as much as possible. However, alternative strategies here, such as relaxed Newton iterations, may prove successful. This is entirely engineering-dependent, dependent on how adaptivity is handled and smoothing that out.
The next stage is when you get to ODEs which are large enough that the sparse factorizations can no longer fit into memory. At this point, you pretty much have to change your nonlinear solver strategy. There are two methods which are being investigated. One is to do Anderson Acceleration. The other is Newton-Krylov methods. Both of these methods cannot just keep something factorized to make proceeding stages in the same step trivial, which means that multi-staged methods (Rosenbrock, ESDIRK, etc.) now have a lot more nonlinear solving time that a single stage method since they have to perform the entire nonlinear solving from scratch every stage. That said, there are one again some engineering factors that get in there. For example, when you implement nonlinear solvers for something like an ESDIRK method, you can use the previous stages to do an extrapolation for the starting value of the nonlinear iterations in the next stage. In many cases, the $c_i$ of the tableau are not monotonic, which in turn means that it's not an extrapolation but a (low order) interpolation, meaning that the nonlinear iteration at some stages converge in one or two goes. So it's not that straightforward to say that multi-staged methods are disadvantaged here without mentioning the accuracy of your stage predictors (which is the term used for this process), so again something "outside the method" and part of the software processes is influencing the practical results. But, in this domain BDF does seem to shine because it only has one stage. However, the fact that it uses such small steps is somewhat of a concern (lack of L-stability above order 2, and the lack of optimized leading truncation error coefficients gives this property) and does lead to quite a few linear solves, so I wouldn't be surprised if well-designed and well-implemented ESDIRKs can topple that sooner or later.
This focused on the methods which have been the most benchmarked, but alternatively the EPRIK and exponential Rosenbrock methods are an entirely alternative branch to explore. There are some good results that have been reported in the literature, and DifferentialEquations.jl does have the full Krylov-expmv + Krylov adaptivity setup in ExponentialUtilities.jl and implemented in 5th order EPIRK methods and so if you want to join our benchmarking team we would love to see how these are turning out. And of course, since engineering is so important, these are quite new and will likely need a few rounds of profiling. Some preliminary results are exciting (what's interesting is that while ARKODE fails on a lot of ODE benchmarks, it's quite good on discretized PDE benchmarks).
What these benchmarks are also highlighting is that one of the ways to improve stiff solvers in 2019 is to not solve it as a single f, but as a split f = f1 + f2 and treat only one part implicitly via IMEX integrators. IMEX BDF (called SBDF), IMEX ESDIRK, IMEX Rosenbrock, etc. all exist and are interesting routes that need more testing.
Also, if you only have "semi-stiffness", then stabilized explicit methods are explicit methods with enough adaptive stability to handle some stiff problems. The implementations exist but need benchmarks. We are currently undergoing a benchmark overhaul to update to use Weave.jl files (for auto-running) and update all benchmarks to Julia v1.0, so stabilized RK, RosW, exponential integrators, and others will be joining the benchmarks once this is done. Our benchmarks also are missing the efficiency when GPUs and TPUs are integrated into the software: the native Julia methods are close to being able to use GPUs so it'll be interesting to test what happens when we make use of it. That puts a higher stress on memory but a lower stress on the cost of linear algebra routines, maybe bumping things more towards the $f$-dominated domain.
Because this is a gigantic topic, there is a chance I also missed one of the infinite options you specifically wanted more discussion on. Oops I'm sorry but it would take $O(n^5)$ time to review all methods. Ask again in 5 years and we'll have more benchmarks :).
So in total, you wanted a short summary about the methods and not the software, but since you asked about what's actually state-of-the-art I had to respond with all of the different ways that the results are actually tied to the software. However, these are a bunch of different directions that seem to be fruitful in different domains, and the JuliaDiffEq team is continuing to implement, optimize, and benchmark. So any time now some of these results may be changing, and that's just how it is since asymptotic factors and pure mathematical constraints only matter at infinity. If you want to know more about the field, summaries like this can be helpful but you will never get a good idea until you get your hands dirty and start benchmarking and profiling. We'd be happy to have you submit problems from your domain to our benchmarks and see what does well!
P.S. To track down the papers on the methods, see our citing page. That's not completely up to date, so you may have to dive into our issues in some cases. A PR to improve our citations page is always welcome of course.
Edit: Yingbo Ma shared this nice little script that highlights the stepping behavior difference of BDF vs the other methods:
using OrdinaryDiffEq, Plots, ParameterizedFunctions, Sundials, ODEInterfaceDiffEq
hires = @ode_def Hires begin
dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
dy2 = 1.71*y1 - 8.75*y2
dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
0.43*y6 + 0.69*y7
dy7 = 280.0*y6*y8 - 1.81*y7
dy8 = -280.0*y6*y8 + 1.81*y7
end
u0 = zeros(8)
u0[1] = 1
u0[8] = 0.0057
prob = ODEProblem(hires,u0,(0.0,321.8122))
function plotdts(sols::Tuple, title; kwargs...)
ts = map(s->s.second.t, sols)
dts = map(t->[0; diff(t)], ts)
xlim = max(map(maximum, ts)...)
plt = plot(title=title, xlabel="t", ylabel="dt", ylims=(0, max(map(maximum, dts)...)*1.3), xlims=(0, xlim*1.5); kwargs...)
vline!(plt, [xlim], lab="t1")
for i in 1:length(ts)
plot!(plt, ts[i], dts[i], lab=sols[i].first, marker=(2, ))
end
plt
end
plotdts(("CVODE_BDF"=>solve(prob, CVODE_BDF()),
"KenCarp4"=>solve(prob, KenCarp4()),
"radau5"=>solve(prob, radau5())), "Time vs Step Size", dpi=200)
savefig("$(homedir())/tvsdt.png")
