Numerical methods for ordinary differential equations are methods used to findnumerical approximations to the solutions ofordinary differential equations (ODEs). Their use is also known as "numerical integration", although this term can also refer to the computation ofintegrals.

Many differential equations cannot be solved exactly. For practical purposes, however – such as in engineering – a numeric approximation to the solution is often sufficient. Thealgorithms studied here can be used to compute such an approximation. An alternative method is to use techniques fromcalculus to obtain aseries expansion of the solution.

Ordinary differential equations occur in many scientific disciplines, includingphysics,chemistry,biology, andeconomics.^[1] In addition, some methods innumerical partial differential equations convert thepartial differential equation into an ordinary differential equation, which must then be solved.

A first-order differential equation is anInitial value problem (IVP) of the form,^[2]

${\displaystyle y^{\prime }(t)=f(t,y(t)),y(t_0)=y_0,}$

1

where${\displaystyle f}$ is a function${\displaystyle f:[t_0,\mathrm{\infty })\times {\mathbb{R}}^d\to {\mathbb{R}}^d}$, and the initial condition${\displaystyle y_0\in {\mathbb{R}}^d}$ is a given vector.First-order means that only the first derivative ofy appears in the equation, and higher derivatives are absent.

Without loss of generality to higher-order systems, we restrict ourselves tofirst-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 equationy′′ = −y can be rewritten as two first-order equations:y′ =z andz′ = −y.

In this section, we describe numerical methods for IVPs, and remark thatboundary value problems (BVPs) require a different set of tools. In a BVP, one defines values, or components of the solutiony at more than one point. Because of this, different methods need to be used to solve BVPs. For example, theshooting method (and its variants) or global methods likefinite differences,^[3]Galerkin methods,^[4] orcollocation methods are appropriate for that class of problems.

ThePicard–Lindelöf theorem states that there is a unique solution, providedf isLipschitz-continuous.

Numerical methods for solving first-order IVPs often fall into one of two large categories:^[5]linear multistep methods, orRunge–Kutta methods. A further division can be realized by dividing methods into those that are explicit and those that are implicit. For example, implicitlinear multistep methods includeAdams-Moulton methods, andbackward differentiation methods (BDF), whereasimplicit Runge–Kutta methods^[6] include diagonally implicit Runge–Kutta (DIRK),^[7][8] singly diagonally implicit Runge–Kutta (SDIRK),^[9] and Gauss–Radau^[10] (based onGaussian quadrature^[11]) numerical methods. Explicit examples from thelinear multistep family include theAdams–Bashforth methods, and any Runge–Kutta method with a lower diagonalButcher tableau isexplicit. A loose rule of thumb dictates thatstiff differential equations require the use of implicit schemes, whereas non-stiff problems can be solved more efficiently with explicit schemes.

The so-calledgeneral linear methods (GLMs) are a generalization of the above two large classes of methods.^[12]

From any point on a curve, you can find an approximation of a nearby point on the curve by moving a short distance along a linetangent to the curve.

Starting with the differential equation (1), we replace the derivativey′ by thefinite difference approximation

${\displaystyle y^{\prime }(t)\approx \frac{y(t+h)-y(t)}{h},}$

2

which when re-arranged yields the following formula

${\displaystyle y(t+h)\approx y(t)+h{y}^{\prime }(t)}$

and using (1) gives:

${\displaystyle y(t+h)\approx y(t)+hf(t,y(t)).}$

3

This formula is usually applied in the following way. We choose a step sizeh, and we construct the sequence${\displaystyle t_0,t_1=t_0+h,t_2=t_0+2h,...}$ We denote by${\displaystyle y_n}$ a numerical estimate of the exact solution${\displaystyle y(t_n)}$. Motivated by (3), we compute these estimates by the followingrecursive scheme

${\displaystyle y_{n+1}=y_n+hf(t_n,y_n).}$

4

This is theEuler method (orforward Euler method, in contrast with thebackward Euler method, to be described below). The method is named afterLeonhard Euler who described it in 1768.

The Euler method is an example of anexplicit method. This means that the new valuey_n+1 is defined in terms of things that are already known, likey_n.

If, instead of (2), we use the approximation

${\displaystyle y^{\prime }(t)\approx \frac{y(t)-y(t-h)}{h},}$

5

we get thebackward Euler method:

${\displaystyle y_{n+1}=y_n+hf(t_{n+1},y_{n+1}).}$

6

The backward Euler method is animplicit method, meaning that we have to solve an equation to findy_n+1. One often usesfixed-point iteration or (some modification of) theNewton–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 (6) is that they are usually more stable for solving astiff equation, meaning that a larger step sizeh can be used.

Exponential integrators describe a large class of integrators that have recently seen a lot of development.^[13] They date back to at least the 1960s.

In place of (1), we assume the differential equation is either of the form

${\displaystyle y^{\prime }(t)=-Ay+\mathcal{N}(y),}$

7

or it has been locally linearized about a background state to produce a linear term${\displaystyle -Ay}$ and a nonlinear term${\displaystyle \mathcal{N}(y)}$.

Exponential integrators are constructed by multiplying (7) by$e^{At}$, and exactly integrating the result overa time interval${\displaystyle [t_n,t_{n+1}=t_n+h]}$:

${\displaystyle y_{n+1}=e^{-Ah}{y}_n+{\int }_0^h{e}^{-(h-\tau )A}\mathcal{N}\left(y\left(t_n+\tau \right)\right)d\tau .}$

This integral equation is exact, but it doesn't define the integral.

The first-order exponential integrator can be realized by holding${\displaystyle \mathcal{N}(y(t_n+\tau ))}$ constant over the full interval:

${\displaystyle y_{n+1}=e^{-Ah}{y}_n+A^{-1}(1-e^{-Ah})\mathcal{N}(y(t_n)).}$

8

The Euler method is often not accurate enough. In more precise terms, it only has order one (the concept oforder is explained below). This caused mathematicians to look for higher-order methods.

One possibility is to use not only the previously computed valuey_n to determiney_n+1, but to make the solution depend on more past values. This yields a so-calledmultistep method. Perhaps the simplest is theleapfrog method which is second order and (roughly speaking) relies on two time values.

Almost all practical multistep methods fall within the family oflinear multistep methods, which have the form

Another possibility is to use more points in the interval${\displaystyle [t_n,t_{n+1}]}$. This leads to the family ofRunge–Kutta methods, named afterCarl Runge andMartin Kutta. One of their fourth-order methods is especially popular.

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, sovariable 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 anerror 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 avariable order method). Methods based onRichardson extrapolation,^[14] such as theBulirsch–Stoer algorithm,^[15][16] 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 pointst₀,t₁,t₂, ...
• event location: finding the times where, say, a particular function vanishes. This typically requires the use of aroot-finding algorithm.
• support forparallel computing.
• when used for integrating with respect to time, time reversibility

Many methods do not fall within the framework discussed here. Some classes of alternative methods are:

• multiderivative methods, which use not only the functionf but also its derivatives. This class includesHermite–Obreschkoff methods andFehlberg methods, as well as methods like theParker–Sochacki method^[17] or Bychkov–Scherbakov method, which compute the coefficients of theTaylor series of the solutiony 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 methods^[18][19] are especially designed for special classes of ODEs (for example,symplectic integrators for the solution ofHamiltonian equations). They take care that the numerical solution respects the underlying structure or geometry of these classes.
• Quantized state systems methods are a family of ODE integration methods based on the idea of state quantization. They are efficient when simulating sparse systems with frequent discontinuities.

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 ofparallel computing.

Early PinT methods (the earliest being proposed in the 1960s)^[20] 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 ofParareal, a flexible, easy-to-use PinT algorithm that is suitable for solving a wide variety of IVPs. The advent ofexascale 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 powerfulsupercomputers. The most popular methods as of 2023 include Parareal, PFASST, ParaDiag, and MGRIT.^[21]

Numerical analysis is not only the design of numerical methods, but also their analysis. Three central concepts in this analysis are:

• convergence: whether the method approximates the solution,
• order: how well it approximates the solution, and
• stability: whether errors are damped out.^[22]

A numerical method is said to beconvergent if the numerical solution approaches the exact solution as the step sizeh goes to 0. More precisely, we require that for every ODE (1) with aLipschitz functionf and everyt^* > 0,

${\displaystyle \underset{h\to 0^{+}}{lim}\underset{n=0,1,\dots ,\lfloor t^{\ast }/h\rfloor }{max}\Vert y_{n,h}-y(t_n)\Vert =0.}$

All the methods mentioned above are convergent.

Suppose the numerical method is

${\displaystyle y_{n+k}=\mathrm{\Psi }(t_{n+k};y_n,y_{n+1},\dots ,y_{n+k-1};h).}$

Thelocal (truncation) error of the method is the error committed by one step of the method. That is, it is the difference between the result given by the method, assuming that no error was made in earlier steps, and the exact solution:

${\displaystyle {\delta }_{n+k}^h=\mathrm{\Psi }\left(t_{n+k};y(t_n),y(t_{n+1}),\dots ,y(t_{n+k-1});h\right)-y(t_{n+k}).}$

The method is said to beconsistent if

${\displaystyle \underset{h\to 0}{lim}\frac{{\delta }_{n+k}^h}{h}=0.}$

The method hasorder${\displaystyle p}$ if

${\displaystyle {\delta }_{n+k}^h=O(h^{p+1})\text{as }h\to 0.}$

Hence a method is consistent if it has an order greater than 0. The (forward) Euler method (4) and the backward Euler method (6) introduced above both have order 1, so they are consistent. Most methods being used in practice attain higher order. Consistency is a necessary condition for convergence^{[citation needed]}, but not sufficient; for a method to be convergent, it must be both consistent andzero-stable.

A related concept is theglobal (truncation) error, the error sustained in all the steps one needs to reach a fixed time${\displaystyle t}$. Explicitly, the global error at time${\displaystyle t}$ is${\displaystyle y_N-y(t)}$ where${\displaystyle N=(t-t_0)/h}$. The global error of a${\displaystyle p}$th order one-step method is${\displaystyle O(h^p)}$; in particular, such a method is convergent. This statement is not necessarily true for multi-step methods.

For some differential equations, application of standard methods—such as the Euler method, explicitRunge–Kutta methods, ormultistep methods (for example, Adams–Bashforth methods)—exhibit instability in the solutions, though other methods may produce stable solutions. This "difficult behaviour" in the equation (which may not necessarily be complex itself) is described asstiffness, and is often caused by the presence of different time scales in the underlying problem.^[23] For example, a collision in a mechanical system like in animpact oscillator typically occurs at much smaller time scale than the time for the motion of objects; this discrepancy makes for very "sharp turns" in the curves of the state parameters.

Stiff problems are ubiquitous inchemical kinetics,control theory,solid mechanics,weather forecasting,biology,plasma physics, andelectronics. One way to overcome stiffness is to extend the notion of differential equation to that ofdifferential inclusion, which allows for and models non-smoothness.^[24][25]

Below is atimeline of some important developments in this field.^[26][27]

• 1768 -Leonhard Euler publishes his method.
• 1824 -Augustin Louis Cauchy proves convergence of the Euler method. In this proof, Cauchy uses the implicit Euler method.
• 1855 - First mention of themultistep methods ofJohn Couch Adams in a letter written byFrancis Bashforth.
• 1895 -Carl Runge publishes the firstRunge–Kutta method.
• 1901 -Martin Kutta describes the popular fourth-orderRunge–Kutta method.
• 1910 -Lewis Fry Richardson announces hisextrapolation method,Richardson extrapolation.
• 1952 -Charles F. Curtiss andJoseph Oakland Hirschfelder coin the termstiff equations.
• 1963 -Germund Dahlquist introducesA-stability of integration methods.

Boundary value problems (BVPs) are usually solved numerically by solving an approximately equivalent matrix problem obtained by discretizing the original BVP.^[28] The most commonly used method for numerically solving BVPs in one dimension is called theFinite Difference Method.^[3] This method takes advantage of linear combinations of point values to constructfinite difference coefficients that describe derivatives of the function. For example, the second-ordercentral difference approximation to the first derivative is given by:

${\displaystyle \frac{u_{i+1}-u_{i-1}}{2h}=u^{\prime }(x_i)+\mathcal{O}(h^2),}$

and the second-ordercentral difference for the second derivative is given by:

${\displaystyle \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}=u^{″}(x_i)+\mathcal{O}(h^2).}$

In both of these formulae,${\displaystyle h=x_i-x_{i-1}}$ is the distance between neighbouringx values on the discretized domain. One then constructs a linear system that can then be solved by standardmatrix methods. For example, suppose the equation to be solved is:

The next step would be to discretize the problem and use linear derivative approximations such as

${\displaystyle u_i^{″}=\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}}$

and solve the resulting system of linear equations. This would lead to equations such as:

${\displaystyle \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}-u_i=0,\mathrm{\forall }i=1,2,3,...,n-1.}$

On first viewing, this system of equations appears to have difficulty associated with the fact that the equation involves no terms that are not multiplied by variables, but in fact this is false. Ati = 1 andn − 1 there is a term involving the boundary values${\displaystyle u(0)=u_0}$ and${\displaystyle u(1)=u_n}$ and since these two values are known, one can simply substitute them into this equation and as a result have a non-homogeneoussystem of linear equations that has non-trivial solutions.

• Courant–Friedrichs–Lewy condition
• Energy drift
• General linear methods
• List of numerical analysis topics#Numerical methods for ordinary differential equations
• Reversible reference system propagation algorithm
• Modelica Language andOpenModelica software

• Bradie, Brian (2006).A Friendly Introduction to Numerical Analysis. Upper Saddle River, New Jersey: Pearson Prentice Hall.ISBN 978-0-13-013054-9.
• J. C. Butcher,Numerical methods for ordinary differential equations,ISBN 0-471-96758-0
• Hairer, E.; Nørsett, S. P.; Wanner, G. (1993).Solving Ordinary Differential Equations. I. Nonstiff Problems. Springer Series in Computational Mathematics. Vol. 8 (2nd ed.). Springer-Verlag, Berlin.ISBN 3-540-56670-8.MR 1227985.
• Ernst Hairer and Gerhard Wanner,Solving ordinary differential equations II: Stiff and differential-algebraic problems, second edition, Springer Verlag, Berlin, 1996.ISBN 3-540-60452-9.
  (This two-volume monograph systematically covers all aspects of the field.)
• Hochbruck, Marlis; Ostermann, Alexander (May 2010). "Exponential integrators".Acta Numerica.19:209–286.Bibcode:2010AcNum..19..209H.CiteSeerX 10.1.1.187.6794.doi:10.1017/S0962492910000048.S2CID 4841957.
• Arieh Iserles,A First Course in the Numerical Analysis of Differential Equations, Cambridge University Press, 1996.ISBN 0-521-55376-8 (hardback),ISBN 0-521-55655-4 (paperback).
  (Textbook, targeting advanced undergraduate and postgraduate students in mathematics, which also discussesnumerical partial differential equations.)
• John Denholm Lambert,Numerical Methods for Ordinary Differential Systems, John Wiley & Sons, Chichester, 1991.ISBN 0-471-92990-5.
  (Textbook, slightly more demanding than the book by Iserles.)

• Joseph W. Rudmin,Application of the Parker–Sochacki Method to Celestial MechanicsArchived 2016-05-16 at the Portuguese Web Archive, 1998.
• Dominique Tournès,L'intégration approchée des équations différentielles ordinaires (1671–1914), thèse de doctorat de l'université Paris 7 - Denis Diderot, juin 1996. Réimp. Villeneuve d'Ascq : Presses universitaires du Septentrion, 1997, 468 p. (Extensive online material on ODE numerical analysis history, for English-language material on the history of ODE numerical analysis, see, for example, the paper books by Chabert and Goldstine quoted by him.)
• Pchelintsev, A.N. (2020). "An accurate numerical method and algorithm for constructing solutions of chaotic systems".Journal of Applied Nonlinear Dynamics.9 (2):207–221.arXiv:2011.10664.doi:10.5890/JAND.2020.06.004.S2CID 225853788.
• kv onGitHub (C++ library with rigorous ODE solvers)
• INTLAB (A library made byMATLAB/GNU Octave which includes rigorous ODE solvers)

• Numerical differential equations
• Ordinary differential equations

• Articles with short description
• Short description matches Wikidata
• All articles with unsourced statements
• Articles with unsourced statements from September 2019
• Webarchive template other archives