Recommendations for ODE solvers for stiff equations

Question

I'm continuing the research of a former Ph.D. student in my group requiring the solution of a system of ODEs. On a technical note, they wrote:

The system of Boltzmann equations behaves numerically stiff. This means that it is advisable to use an implicit method for its numerical solution to achieve acceptable step sizes (and hence acceptable execution times and numerical errors). Here CVODE with its backward differentiation formula with Newton iteration was used as ODE solver. The full Jacobian was computed analytically in every external step.

I don't have access to their code and I'm new in C/C++, so I wonder if

1. Could you recommend a simpler ODE solver to start with?
2. Is there a tutorial/book/example to learn how to use CVODE? I have looked at the official documentation but I'm afraid it is too high-level for me yet.

More details can be provided if needed. Thanks!

Answer 1

This is a huge open ended question, but I'll copy the recommendation section from the current release of DifferentialEquations.jl:

Stiff Problems

For stiff problems at high tolerances (>1e-2?) it is recommended that you use Rosenbrock23 or TRBDF2. These are robust to oscillations and massive stiffness, though are only efficient when low accuracy is needed. Rosenbrock23 is more efficient for small systems where re-evaluating and re-factorizing the Jacobian is not too costly, and for sufficiently large systems TRBDF2 will be more efficient.

At medium tolerances (>1e-8?) it is recommended you use Rodas5, Rodas4 (the former is more efficient but the later is more reliable), Kvaerno5, or KenCarp4. As native DifferentialEquations.jl solvers, many Julia numeric types (such as BigFloats, ArbFloats, or DecFP) will work.

For faster solving at low tolerances (<1e-9) but when Vector{Float64} is used, use radau.

For asymptotically large systems of ODEs (N>1000?) where f is very costly and the complex eigenvalues are minimal (low oscillations), in that case QNDF will be the most efficient but requires Vector{Float64}. QNDF will also do surprisingly well if the solution is smooth. However, this method can handle less stiffness than other methods and its Newton iterations may fail at low accuracy situations. Other choices to consider in this regime are CVODE_BDF and lsoda.

Special Properties of Stiff Integrators

ImplicitMidpoint is a symmetric and symplectic integrator. Trapezoid is a symmetric (almost symplectic) integrator with adaptive timestepping. ImplicitEuler is an extension to the common algorithm with adaptive timestepping and efficient quasi-Newton Jacobian re-usage which is fully strong-stability preserving (SSP) for hyperbolic PDEs.

Notice that Rodas4 loses accuracy on discretizations of nonlinear parabolic PDEs, and thus it's suggested you replace it with Rodas4P in those situations which is 3rd order. ROS3P is only third order and achieves 3rd order on such problems and can thus be more efficient in this case.

(specifically it's this page)

These recommendations are all supported by the SciMLBenchmarks which are a continuously updating cross-language benchmarking environment for scientific machine learning activities, including differential equation solving and parameter estimation.

Additionally, there are other properties to specialize on. If the system is of the form u'=A(t)u for example, if can use Magnus methods which will be cheap and stable. See this page for a list of methods on specialized forms of non-autonomous linear ODEs. Similarly, there are cases where only some of the equation is stiff, or part of the equation is stiff and linear, in which case split methods like IMEX or exponential integrators are a good idea. See this page for descriptions of IMEX or exponential integrators.

In short, the summary is:

1. Rosenbrock methods are much more efficient for smaller systems but require a more accurate Jacobian. You really only want to be using them if you are using automatic differentiation or symbolic differentiation is used for the Jacobian (you can get away with more lenience if it's a W-moethod, but not much), and if your system is sufficiently small. If those cases hold though, they are by far the fastest.

2. L-stable methods are generally the most efficient, with higher order SDIRK and adaptive time adaptive order BDF methods being good choices. The pure Julia QNDF tends to outperform CVODE_BDF in many cases (related to some details and linear solvers) for the BDF part, and trying TRBDF2 or KenCarp methods is usually a good idea.

3. Use higher order methods when you want lower tolerance. The highest order methods for stiff ODEs are the adaptive order Radau methods.

---

Using Stiff ODE solvers

This is almost a different question from the recommendation, but there are many tutorials in the DifferentialEquaitons.jl documentation to get started by copy and pasting. There is an entire tutorial dedicated to stiff ODEs. So to end, this will solve the Robertson equation with CVODE_BDF:

using DifferentialEquations
function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  k₂*y₂^2
  nothing
end
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4])
sol = solve(prob,CVODE_BDF()) # Don't pass CVODE_BDF() and it will auto-select a method
plot(sol,tspan=(1e-2,1e5),xscale=:log10)

That is probably the easiest way to CVODE up and running as the installation process automates the build of all of the associated libraries and includes the building of LAPACK/BLAS.

Answer 2

CVODE is part of the SUNDIALS package, which is open source. It comes with a manual that has several hundred pages and many example codes. It is excellent software, you should use it!