How can I choose the best algorithm to integrate ODEs numerically?

Question

I have studied in a course several algorithms to integrate ODEs numerical: Runge-Kutta, Predictor-Corrector methods, Taylor...

However the teacher failed to show which is the best for every particular situation. The only thing I know is that implicit methods are appropriate for stiff systems. But how do they compare Runge-Kutta (of any order) with predictor-corrector methods or with Taylor method? Which is best for each situation?

Heuristic answers based on experience may be good enough!

Answer 1

There are many factors to take into account when choosing a numerical integrator. Overall, there are two types of goals for such solutions:

1. Obtain a solution with sufficiently small error in a sufficiently short time (or with limited available computational resources).
2. Ensure that the solution satisfies important qualitative constraints. These may include things like positivity (for densities, energies, and such) conservation (of mass, momentum, energy, or other physically conserved quantities) dissipativity or contractivity.

It is extremely difficult to give general recommendations, but here are some fairly well-agreed-upon suggestions:

• The greater the required accuracy, the more higher order methods become efficient.
• In the end, stiffness depends on the relation between the desired accuracy and the Lipschitz constant for the problem. For low accuracy (or large Lipschitz constant), implicit methods are more efficient. Depending on the spectrum of your problem, you may need L-stability, A-stability, or just $A(0)$-stability.
• For problems where conservation is important: any method is conservative up to discretization errors. If you need conservation to within roundoff error, use a symplectic method.
• If dissipativity or contractivity in an inner-product norm is important, use a B-stable method. If you need one of these properties with respect to a norm that is not induced by an inner product (e.g. $L_1$ or $L_\infty$), strong stability preserving (contractive) methods may be useful.
• If memory is the limiting factor and you need high accuracy (so you wish to use a high-order method), then low-storage Runge-Kutta methods may be useful.

For many problems, the constraints above may still indicate that reasonable choices exist within the class of Runge-Kutta methods or linear multistep methods. At that point the difference becomes an engineering question and may depend on implementation details. You can find useful comparisons between the two for instance in the books of Hairer & Wanner.

There are always exceptions, but...if you have a small non-stiff system (the definition of small depends on your computer and whether the Jacobian is dense or sparse), then usually any explicit integrator of order greater than one will do. If you have a small stiff system, usually any implicit integrator of order greater than one will do. You only need to sweat it if you have a very large system (coming from a PDE or many-body problem) or your ODEs have some unusual feature.

Taylor methods are often reasonable to implement for easy problems, but there is no compelling reason to do so. For hard (i.e. large) problems, Taylor methods can be difficult or impossible to implement, though automatic differentiation may eventually change that.