How can I choose the best algorithm to integrate ODE's numerically?

I have studied in a course several algorithms to integrate ODE's 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!

• What makes you even think that there is an answer to such an question? It's ridiculous to think that there should even be a classification of "every particular situation" let alone what the best method in every case should be. – Johannes Hahn Jun 10 '14 at 23:39
• I don't see it so ridiculous. If there is not an answer then, why doesn't everybody use the same method? For example Runge-Kuttas of different order? What is the point of programming Taylor which is more difficult to program? – Ambesh Jun 11 '14 at 8:49
• Indeed, it is not a ridiculous but an excellent question. However, Johannes Hahn has a point in saying that it is highly unlikely that you will get a satisfactory answer. I've been in numerical methods for more than 30 years and the only thing I've seen is that people put all those different methods in a software package so that other people can make their choice. There seems to be no common consensus about what method is the best, not even when all boundary conditions are known. Maybe saying that Numerical Analysis is just a mess is not even far besides the truth. – Han de Bruijn Jun 11 '14 at 9:31
• @HandeBruijn Well, let's hope for some luck. I accept any partial answer even if based only on experience. – Ambesh Jun 11 '14 at 9:56
• by numerical method, do u mean by using power series – Holy cow Jun 12 '14 at 14:02

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.