Many differential equations cannot be solved analytically, in which case we have to satisfy ourselves with an approximation to the solution. 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 mechanics, chemistry, ecology, and economics. In addition, some methods in numerical partial differential equations convert the partial differential equation in an ordinary differential equation, which must then be solved.
| Table of contents |
|
2 Methods
2.1 The Euler method
3 Analysis2.2 The backward Euler method 2.3 Generalisations 2.4 Advanced features 2.5 Alternative methods 4 History 5 References |
We want to approximate the solution of the differential equation
The above formulation is called an initial value problem
(IVP). The Picard-Lindelöf theorem states that there is a
unique solution, if f is Lipschitz continuous. In contrast,
boundary value problems (BVPs) specify (components of) the
solution y at more than one points. Different methods need to be
used to solve BVPs, for example the shooting method,
multiple shooting or global methods like finite differences or
collocation.
Note that we restrict ourselves to first-order differential
equations (meaning that only the first derivative of y appears in
the equation, and no higher derivatives). However, a higher-order
equation can easily be converted to a first-order equation 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.
Two elementary methods are discussed to give the reader a feeling for
the subject. After that, pointers are provided to other methods (which
are generally more accurate and efficient). The methods mentioned here
are analysed in the next section.
Starting with the differential equation (1), we replace the derivative
y' by the finite difference approximation
The problem
where f is a function that maps
[t0,∞) × Rd
to Rd, and the initial condition
y0 ∈ Rd is a given
vector.Methods
The Euler method
which yields the following formula
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
This is the Euler method, named after Leonhard Euler who
described this method in 1768.
If, instead of (2), we use the approximation
The backward Euler method
we get the backward Euler method:
The backward Euler method is an implicit method, meaning than we
have to solve an equation to find yn+1. One often
uses functional iteration or (some modification of) the
Newton-Raphson method to achieve this. Of course,
it costs time to solve this equation; this cost must be taken into
consideration when one selects the method to use.
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. Almost all practical multistep methods fall within
the family of linear multistep methods, which have the form
Generalisations
Another possibility is to use more points in the interval
[tn,tn+1]. This leads to the
family of Runge-Kutta methods, named after Carle Runge and
Martin Kutta. One of their fourth-order methods is especially
popular.
Both ideas can also be combined. The resulting methods are called general linear methods.
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). Extrapolation methods are often used
to construct various methods of different orders.
Other desirable features include:
Many methods do not fall within the framework discussed here. Some
classes of alternative methods are:
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).
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,
Suppose the numerical method is
The local error is the error committed in a single step. A related
concept is the global error, the error sustained in all the steps
one needs to reach a fixed time t. Explicitly, the global error at
time t is
y(t-t0)/h - y(t).
The global error of a pth order one-step method (that is, a method
of the form (4) with k = 1) is O(hp);
in particular, such a method is convergent. This statement is not
necessarily true for multi-step methods.
Loosely speaking, a numerical method is called
stable if unwanted components in the numerical
solution die out over time. Many different aspects of stability have
been discussed in the literature. We will only treat one of them.
A method is A-stable if the numerical results yn
approach zero as n → 0 for all values of the step
size h when this method is applied to the equation
y' = λy for all
λ ∈ C with
Re λ < 0. Note that for this equation, the exact
solution also goes to zero. The (forward) Euler method is not
A-stable, but the backward Euler method is A-stable.
For some differential equations, it does not matter much whether the
method is stable. However, for other equations, stable methods perform
far better; these equations are said to be stiff (it is hard to
formulate a more precise definition). Stiffness is often caused by the
presence of different time scales in the underlying problem. Stiff
problems are ubiquitous in (chemical) kinetics, control theory,
weather prediction, biology, and electronics.
Below is a concise timeline of some important developments in this
field.
Advanced features
Alternative methods
Analysis
Convergence
All the methods mentioned above are convergent. In fact, convergence
is a condition sine qua non for any numerical scheme.Order
The method is said to have order p if
The quantity on the left-hand side is called the local error of
the method. The (forward) Euler method and the backward Euler method
introduced above both have order 1. Most methods being used in
practise attain higher order.Stability and stiffness
History
References