Numerical methods for ordinary differential equations

Illustration of numerical integration for the differential equation y'=y, y(0)=1. Blue: the Euler method, green: the midpoint method, red: the exact solution, y=e^t. The step size is h=1.0.
The same illustration for h=0.25. It is seen that the midpoint method converges faster than the Euler method.

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 is sometimes taken to mean the computation of integrals.

Many differential equations cannot be solved using symbolic computation ("analysis"). For practical purposes, however – such as in engineering – a numeric 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, for instance in physics, chemistry, biology, and economics. In addition, some methods in numerical partial differential equations convert the partial differential equation into an ordinary differential equation, which must then be solved.

The problem

A first-order differential equation is an Initial value problem (IVP) of the form,[1]

y'(t) = f(t,y(t)), \qquad y(t_0)=y_0, \qquad\qquad (1)

where f is a function that maps [t0,∞) × Rd to Rd, and the initial condition y0  Rd is a given vector. First-order means that only the first derivative of y 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 y'' = y can be rewritten as two first-order equations: y' = z and z' = y.

In this section, we describe numerical methods for IVPs, and remark that boundary value problems (BVPs) require a different set of tools. In a BVP, one defines values, or components of the solution y at more than one point. Because of this, different methods need to be used to solve BVPs. For example, the shooting method (and its variants) or global methods like finite differences, Galerkin methods, or collocation methods are appropriate for that class of problems.

The Picard–Lindelöf theorem states that there is a unique solution, provided f is Lipschitz-continuous.

Methods

Numerical methods for solving first-order IVPs often fall into one of two large categories: linear multistep methods, or Runge-Kutta methods. A further division can be realized by dividing methods into those that are explicit and those that are implicit. For example, implicit linear multistep methods include Adams-Moulton methods, and backward differentiation methods (BDF), whereas implicit Runge-Kutta methods[2] include diagonally implicit Runge-Kutta (DIRK), singly diagonally implicit runge kutta (SDIRK), and Gauss-Radau (based on Gaussian quadrature) 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.

Euler method

For more details on this topic, see Euler method.

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 line tangent to the curve.

Starting with the differential equation (1), we replace the derivative y' by the finite difference approximation

 y'(t) \approx \frac{y(t+h) - y(t)}{h}, \qquad\qquad (2)

which when re-arranged yields the following formula

 y(t+h) \approx y(t) + hy'(t) \qquad\qquad

and using (1) gives:

 y(t+h) \approx y(t) + hf(t,y(t)). \qquad\qquad (3)

This formula is usually applied in the following way. We choose a step size h, and we construct the sequence t0, t1 = t0 + h, t2 = t0 + 2h, … We denote by yn a numerical estimate of the exact solution y(tn). Motivated by (3), we compute these estimates by the following recursive scheme

 y_{n+1} = y_n + hf(t_n,y_n). \qquad\qquad (4)

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 yn+1 is defined in terms of things that are already known, like yn.

Backward Euler method

For more details on this topic, see Backward Euler method.

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

 y'(t) \approx \frac{y(t) - y(t-h)}{h}, \qquad\qquad (5)

we get the backward Euler method:

 y_{n+1} = y_n + hf(t_{n+1},y_{n+1}). \qquad\qquad (6)

The backward Euler method is an implicit method, meaning that we have to solve an equation to find yn+1. 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 (6) is that they are usually more stable for solving a stiff equation, meaning that a larger step size h can be used.

First-order exponential integrator method

For more details on this topic, see Exponential integrators.

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

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

y'(t) = -A\, y+ \mathcal{N}(y), \qquad\qquad\qquad (7)

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

Exponential integrators are constructed by multiplying (7) by e^{A t}, and exactly integrating the result over a time interval [t_n, t_{n+1} = t_n + h]:

 y_{n+1} = e^{-A h } y_n + \int_{0}^{h} e^{ -(h-\tau) A } \mathcal{N}\left( y\left( t_n+\tau \right) \right)\, d\tau.

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

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

y_{n+1} = e^{-Ah}y_n + A^{-1}(1-e^{-Ah}) \mathcal{N}( y( t_n ) )\ . \qquad\qquad (8)

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

 \alpha_k y_{n+k} + \alpha_{k-1} y_{n+k-1} + \cdots 
+ \alpha_0 y_n
 = h \left[ \beta_k f(t_{n+k},y_{n+k}) + \beta_{k-1}
f(t_{n+k-1},y_{n+k-1}) + \cdots + \beta_0 f(t_n,y_n) \right].

Another possibility is to use more points in the interval [tn,tn+1]. 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, such as the Bulirsch–Stoer algorithm, are often used to construct various methods of different orders.

Other desirable features include:

Alternative methods

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

Parallel-in-time methods

For applications that require parallel computing on supercomputers, the degree of concurrency offered by a numerical method becomes relevant. In view of the challenges from exascale computing systems, numerical methods for initial value problems which can provide concurrency in temporal direction are being studied.[4] Parareal is a relatively well known example of such a parallel-in-time integration method, but early ideas go back into the 1960s.[5]

Analysis

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

Convergence

A numerical method is said to be convergent if the numerical solution approaches the exact solution as the step size h goes to 0. More precisely, we require that for every ODE (1) with a Lipschitz function f and every t* > 0,

 \lim_{h\to0+} \max_{n=0,1,\dots,\lfloor t^*/h\rfloor} \| y_{n,h} - y(t_n) \| = 0.

All the methods mentioned above are convergent. In fact, a numerical scheme has to be convergent to be of any use.

Consistency and order

For more details on this topic, see Truncation error (numerical integration).

Suppose the numerical method is

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

The local (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:

 \delta^h_{n+k} = \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 be consistent if

 \lim_{h\to 0} \frac{\delta^h_{n+k}}{h} = 0.

The method has order p if

 \delta^h_{n+k} = O(h^{p+1}) \quad\mbox{as } h\to0.

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, but not sufficient; for a method to be convergent, it must be both consistent and zero-stable.

A related concept is the global (truncation) error, the error sustained in all the steps one needs to reach a fixed time t. Explicitly, the global error at time t is yN  y(t) where N = (tt0)/h. The global error of a pth order one-step method is O(hp); in particular, such a method is convergent. This statement is not necessarily true for multi-step methods.

Stability and stiffness

For more details on this topic, see Stiff equation.

For some differential equations, application of standard methods such as the Euler method, explicit Runge–Kutta methods, or multistep methods (e.g., 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 as stiffness, and is often caused by the presence of different time scales in the underlying problem. For example, a collision in a mechanical system like in an impact 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 in chemical kinetics, control theory, solid mechanics, weather forecasting, biology, plasma physics, and electronics. One way to overcome stiffness is to extend the notion of differential equation to that of differential inclusion, which allows for and models non-smoothness.[6][7]

History

Below is a timeline of some important developments in this field.

Numerical solutions to second-order one-dimensional boundary value problems

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

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

and the second-order central difference for the second derivative is given by:

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

In both of these formulae,  h=x_i-x_{i-1} is the distance between neighbouring x values on the discretized domain. One then constructs a linear system that can then be solved by standard matrix methods. For instance, suppose the equation to be solved is:

 \frac{d^2 u}{dx^2} -u =0,
 u(0)=0,
 u(1)=1.

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

 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:

 \frac{u_{i+1}-2u_{i}+u_{i-1}}{h^2}-u_i = 0, \quad \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. At i = 1 and n 1 there is a term involving the boundary values u(0)=u_0 and  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-homogeneous linear system of equations that has non-trivial solutions.

See also

Notes

  1. Bradie (2006, pp. 533–655)
  2. Hairer, Nørsett & Wanner (1993), pp. 204–215
  3. Hochbruck (2010, pp. 209–286) This is a modern and extensive review paper for exponential integrators
  4. Gander, Martin J. (2015). 50 years of Time Parallel Time Integration. Contributions in Mathematical and Computational Sciences (1 ed.). Springer International Publishing. doi:10.1007/978-3-319-23321-5. ISBN 978-3-319-23321-5.
  5. Nievergelt, Jürg (1964). "Parallel methods for integrating ordinary differential equations". Communications of the ACM (ACM New York) 7 (12): 731–733. doi:10.1145/355588.365137. Retrieved August 2015.
  6. Markus Kunze and Tassilo Kupper (2001). "Non-smooth Dynamical Systems: An Overview". In Bernold Fiedler. Ergodic Theory, Analysis, and Efficient Simulation of Dynamical Systems. Springer Science & Business Media. p. 431. ISBN 978-3-540-41290-8.
  7. Thao Dang (2011). "Model-Based Testing of Hybrid Systems". In Justyna Zander, Ina Schieferdecker and Pieter J. Mosterman. Model-Based Testing for Embedded Systems. CRC Press. p. 411. ISBN 978-1-4398-1845-9.

References

External links

This article is issued from Wikipedia - version of the Friday, March 04, 2016. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.