# What is the state of the art in solving stiff initial value problems?

I'm looking for current references on solving stiff ODEs. Most of what I know (say, BDF methods) apparently date back to the 1980's, and I feel like a lot of progress should have been made in that time, but I'm having trouble finding the "right" search terms to get what I need. What are the (current) best methods for solving stiff IVPs?

(Note: I'm not looking for software, but rather the methods themselves. I've seen a lot of people around the web recommending to "just use MATLAB", but I need something for a very specific application that MATLAB would not be able to handle anyway.)

• There is no singular state of the art for solving stiff ODEs; the field of study is much too wide. If you can specify your needs a bit, chances are you will get more useful answers. What kind of application do you have in mind? What is it about the problem that prevents the use of MATLAB? Are you looking for fully implicit methods or do you need something that might be easier to implement / faster / different for some other reason? Is accuracy / speed / simplicity / robustness / something else a priority? Does the problem have some regularity that should be preserved by the solver? – ekkilop Apr 6 at 9:46
• Julia has implemented more modern methods, I would say that their differential equations package is the closest to the state of the art that you would get. – nicoguaro Apr 6 at 12:04
• Something Matlab couldn't handle? I doubt there exists such a thing. – mathreadler Apr 6 at 21:42

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()),
savefig("$(homedir())/tvsdt.png")  • Whoa, this is an amazing answer! I'll give a bit more detail about my aversion to MATLAB. I am using vectors and matrices defined in a custom C++ format, so there's the issue of language, but also I am having trouble finding any library that will let me define a custom vector implementation. CVODE does do this, but some of the operations it expects are difficult (impossible?) to implement on this vector space (I'm thinking specifically of vector element-wise division, for example). I thought perhaps I could write a general solver given a method for very stiff, chemical kinetics equations. – emprice Apr 6 at 13:21 • Vector element-wise division likely only occurs in CVODE's error estimation phase, and they do allow you to override the error estimator. Making that something simple like an L1 difference might make that a bit easier (it would also get rid of the difficulty of defining element-wise square root, which would come next :)). That might be all you need to get a working N_Vector for your problem, if your N_Vector is compatible with the linear algebra routines you require. – Chris Rackauckas Apr 6 at 13:25 • Yes, "general solver given a method for very stiff, chemical kinetics equations" doesn't really exist and depends on how stiff, how many ODEs, cost of the chemical reaction equations (how many nonlinear feedbacks do you have? Those are usually the nontrivial terms), etc. all matter for what method ends up being efficient in the end. So when in doubt, when using C++, throwing Sundials at the problem is usually a good bet. It takes quite a bit of engineering to overcome LSODX/Sundials/DifferentialEquations.jl/Hairer on most problems these days, given the fine tuning of jac-reuse and stepping. – Chris Rackauckas Apr 6 at 13:28 • Thanks for the info! I didn't know it could be something that simple -- perhaps I should have just asked about N_Vector implementations in the first place. I did think it was odd to require that operation since it's so specialized and only really makes sense if you're using a more traditional vector space. – emprice Apr 6 at 13:31 • BDF itself only requires doing linear combinations (if you don't divide by$gamma\$ in the W implementation, which I don't think Sundials does). But the error estimator calculates the difference of two steps as E and then relates that to the relative and absolute tolerances, the relative tolerances making E2 sized the same as the system since it uses the values of the current state, and then checks if ||E./E2|| < 1, basically asking whether the error is lower than the tolerance. So unless there's some other funny business going on, that's the only place it should occur. – Chris Rackauckas Apr 6 at 13:34