Numerical methods for ordinary differential equations are methods used to find numerical approximations to the solutions of ordinary differential equations (ODEs). Their use is also known as "numerical integration", although this term can also refer to the computation of integrals.
Many differential equations cannot be solved exactly. For practical purposes, however, a numerical approximation to the solution is often sufficient. The algorithms studied here can be used to compute such an approximation. An alternative method is to use techniques from calculus to obtain a series expansion of the solution.
Ordinary differential equations occur in many scientific disciplines, including physics, chemistry, biology, and economics.Chicone, C. (2006). Ordinary differential equations with applications (Vol. 34). Springer Science & Business Media. In addition, some numerical methods for partial differential equations convert the partial differential equation into a system of ordinary differential equations, which can then be solved numerically.
Initial value problems
An Initial value problem (IVP) is a first-order differential equation plus an initial condition,
where is a function , and is a given vector. First-order means that only the first derivative of appears in the equation, and higher derivatives are absent.
Without loss of generality to higher-order systems, we restrict ourselves to first-order differential equations, because a higher-order ODE can be converted into a larger system of first-order equations by introducing extra variables. For example, the second-order equation can be rewritten as two first-order equations: and .
In this section, we describe numerical methods for IVPs. Boundary value problems (BVPs) require a different set of tools, discussed in the section below.
The Picard–Lindelöf theorem states that the IVP has a unique solution, provided that is Lipschitz-continuous.
Methods
Numerical methods for solving first-order IVPs often fall into one of two large categories:Griffiths, D. F., & Higham, D. J. (2010). Numerical methods for ordinary differential equations: initial value problems. Springer Science & Business Media. linear multistep methods, which make use of information from previous steps to obtain the next step, or Runge–Kutta methods, which use intermediate stages to calculate the next step.
A further division is between explicit and implicit methods. Explicit methods only involve information from the current and previous steps, whereas implicit methods involve the next step, requiring the solution of an equation to solve for the next step.
For example, implicit linear multistep methods include Adams-Moulton methods, and backward differentiation methods (BDF), whereas implicit Runge–Kutta methods include diagonally implicit Runge–Kutta (DIRK),Alexander, R. (1977). Diagonally implicit Runge–Kutta methods for stiff ODE's. SIAM Journal on Numerical Analysis, 14(6), 1006-1021.Cash, J. R. (1979). Diagonally implicit Runge-Kutta formulae with error estimates. IMA Journal of Applied Mathematics, 24(3), 293-301. singly diagonally implicit Runge–Kutta (SDIRK),Ferracina, L., & Spijker, M. N. (2008). Strong stability of singly-diagonally-implicit Runge–Kutta methods. Applied Numerical Mathematics, 58(11), 1675-1686. and Gauss–RadauEverhart, E. (1985). An efficient integrator that uses Gauss-Radau spacings. In International Astronomical Union Colloquium (Vol. 83, pp. 185–202). Cambridge University Press. (based on Gaussian quadratureWeisstein, Eric W. "Gaussian Quadrature." From MathWorld--A Wolfram Web Resource. https://mathworld.wolfram.com/GaussianQuadrature.html) numerical methods. Explicit examples from the linear multistep family include the Adams–Bashforth methods, and any Runge–Kutta method with a lower diagonal Butcher tableau is explicit. A loose rule of thumb dictates that stiff differential equations require the use of implicit schemes, whereas non-stiff problems can be solved more efficiently with explicit schemes.
The so-called general linear methods (GLMs) are a generalization of the above two large classes of methods.Butcher, J. C. (1987). The numerical analysis of ordinary differential equations: Runge-Kutta and general linear methods. Wiley-Interscience.
Euler method
Euler's method can be thought of as approximating a nearby point on a curve by moving a short distance along a line tangent to the curve.
Starting with the differential equation (), we replace the derivative by the finite difference approximation
which when re-arranged yields the following formula
and using () gives:
This formula is usually applied in the following way. We choose a step size h, and we construct the sequence We denote by a numerical estimate of the exact solution . Motivated by (), we compute these estimates by the following recursive scheme
This is the Euler method (or forward Euler method, in contrast with the backward Euler method, to be described below). The method is named after Leonhard Euler who described it in 1768.
The Euler method is an example of an explicit method. This means that the new value is defined in terms of things that are already known, like .
Backward Euler method
If, instead of (), we use the approximation
we get the backward Euler method:
The backward Euler method is an implicit method, meaning that we have to solve an equation to find . One often uses fixed-point iteration or (some modification of) the Newton–Raphson method to achieve this.
It costs more time to solve this equation than explicit methods; this cost must be taken into consideration when one selects the method to use. The advantage of implicit methods such as () is that they are usually more stable for solving a stiff equation, meaning that a larger step size can be used.
First-order exponential integrator method
Exponential integrators describe a large class of integrators that have recently seen a lot of development. This is a modern and extensive review paper for exponential integrators They date back to at least the 1960s.
In place of (), we assume the differential equation is either of the form
or it has been locally linearized about a background state to produce a linear term and a nonlinear term .
Exponential integrators are constructed by multiplying () by , and exactly integrating the result over a time interval where
This integral equation is exact, but it does not define the integral.
The first-order exponential integrator can be obtained by approximating as a constant over the full interval, so that the integral can be evaluated:
Exponential integrators are particularly useful for stiff differential equations where the stiffness lies in the linear term. The linear term is solved exactly, so it does not limit the step size. Higher order exponential integrators can be derived, either of multistep type or Runge-Kutta type.
Generalizations The Euler method is often not accurate enough. In more precise terms, it only has order one (the concept of order is explained below). This caused mathematicians to look for higher-order methods.
One possibility is to use not only the previously computed value yn to determine yn+1, but to make the solution depend on more past values. This yields a so-called multistep method. Perhaps the simplest is the leapfrog method which is second order and (roughly speaking) relies on two time values.
Almost all practical multistep methods fall within the family of linear multistep methods, which have the form
Another possibility is to use more points in the interval . This leads to the family of Runge–Kutta methods, named after Carl Runge and Martin Kutta. One of their fourth-order methods is especially popular.
Advanced features A good implementation of one of these methods for solving an ODE entails more than the time-stepping formula.
It is often inefficient to use the same step size all the time, so variable step-size methods have been developed. Usually, the step size is chosen such that the (local) error per step is below some tolerance level. This means that the methods must also compute an error indicator, an estimate of the local error.
An extension of this idea is to choose dynamically between different methods of different orders (this is called a variable order method). Methods based on Richardson extrapolation,Brezinski, C., & Zaglia, M. R. (2013). Extrapolation methods: theory and practice. Elsevier. such as the Bulirsch–Stoer algorithm,Monroe, J. L. (2002). Extrapolation and the Bulirsch-Stoer algorithm. Physical Review E, 65(6), 066116.Kirpekar, S. (2003). Implementation of the Bulirsch Stoer extrapolation method. Department of Mechanical Engineering, UC Berkeley/California. are often used to construct various methods of different orders.
Other desirable features include: dense output: cheap numerical approximations for the whole integration interval, and not only at the points t0, t1, t2, ... event location: finding the times where, say, a particular function vanishes. This typically requires the use of a root-finding algorithm. support for parallel computing. when used for integrating with respect to time, time reversibility
Alternative methods Many methods do not fall within the framework discussed here. Some classes of alternative methods are: multiderivative methods, which use not only the function f but also its derivatives. This class includes Hermite–Obreschkoff methods and Fehlberg methods, as well as methods like the Parker–Sochacki methodNurminskii, E. A., & Buryi, A. A. (2011). Parker-Sochacki method for solving systems of ordinary differential equations using graphics processors. Numerical Analysis and Applications, 4(3), 223. or Bychkov–Scherbakov method, which compute the coefficients of the Taylor series of the solution y recursively. methods for second order ODEs. We said that all higher-order ODEs can be transformed to first-order ODEs of the form (1). While this is certainly true, it may not be the best way to proceed. In particular, Nyström methods work directly with second-order equations. geometric integration methodsHairer, E., Lubich, C., & Wanner, G. (2006). Geometric numerical integration: structure-preserving algorithms for ordinary differential equations (Vol. 31). Springer Science & Business Media.Hairer, E., Lubich, C., & Wanner, G. (2003). Geometric numerical integration illustrated by the Störmer–Verlet method. Acta Numerica, 12, 399-450. are especially designed for special classes of ODEs (for example, symplectic integrators for the solution of Hamiltonian equations). They take care that the numerical solution respects the underlying structure or geometry of these classes.
Parallel-in-time methods Some IVPs require integration at such high temporal resolution and/or over such long time intervals that classical serial time-stepping methods become computationally infeasible to run in real-time (e.g. IVPs in numerical weather prediction, plasma modelling, and molecular dynamics). Parallel-in-time (PinT) methods have been developed in response to these issues in order to reduce simulation runtimes through the use of parallel computing.
Early PinT methods (the earliest being proposed in the 1960s) were initially overlooked by researchers due to the fact that the parallel computing architectures that they required were not yet widely available. With more computing power available, interest was renewed in the early 2000s with the development of Parareal, a flexible, easy-to-use PinT algorithm that is suitable for solving a wide variety of IVPs. The advent of exascale computing has meant that PinT algorithms are attracting increasing research attention and are being developed in such a way that they can harness the world's most powerful supercomputers. The most popular methods as of 2023 include Parareal, PFASST, ParaDiag, and MGRIT.
Analysis