Numerical integration

In numerical analysis, numerical integration constitutes a broad family of algorithms for calculating the numerical value of a definite integral, and by extension, the term is also sometimes used to describe numerical algorithms for solving differential equations.

Methods for one-dimensional integrals

Numerical integration methods can generally be described as combining evaluations of the integrand to get an approximation to the integral.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

An important part of the analysis of any numerical integration method is to study the behavior of the approximation error as a function of the number of integrand evaluations.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

A method which yields a small error for a small number of evaluations is usually considered superior.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Reducing the number of evaluations of the integrand reduces the number of arithmetic operations involved,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

and therefore reduces the total round-off error.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Also,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

each evaluation takes time, and the integrand may be arbitrarily complicated.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

It should be noted, however, that a 'brute force' kind of numerical integration can always be done, in a very simplistic way, by evaluating the integrand with very small increments.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Quadrature rules based on interpolating functions

A large class of quadrature rules can be derived by constructing interpolating functions which are easy to integrate.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Typically these interpolating functions are polynomials.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

The simplest method of this type is to let the interpolating function be a constant function (a polynomial of order zero)

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

which passes through the point ((a+b)/2, f((a+b)/2)).

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

This is called the midpoint rule or rectangle rule.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

:int_a^b f(x),dx pprox (b-a) , fleft(rac{a+b}{2} ight)

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

The interpolating function may be an affine function (a polynomial of degree 1)

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

which passes through the points (a, f(a)) and (b, f(b)).

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

This is called the trapezoidal rule.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

:int_a^b f(x),dx pprox (b-a) , rac{f(a) + f(b)}{2}

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

For either one of these rules,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

we can make a more accurate approximation by breaking up the interval into some number n of subintervals,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

computing an approximation for each subinterval,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

then adding up all the results.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

This is called a composite rule, extended rule, or iterated rule.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

For example,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

the composite trapezoidal rule can be stated as

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

:int_a^b f(x),dx pprox rac{b-a}{n} left( {f(a) + f(b) over 2} + sum_{k=1}^{n-1} f left( a+k rac{b-a}{n} ight) ight)

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

where the subintervals have the form , with h = (b-a)/n and k = 0, 1, 2, ..., n-1.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Interpolation with polynomials evaluated at equally-spaced points in yields the Newton-Cotes formulas,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

of which the rectangle rule and the trapezoidal rule are examples.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Simpson's rule, which is based on a polynomial of order 2, is also a Newton-Cotes formula.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

If we allow the intervals between interpolation points to vary, we find another group of quadrature formulas, called Gaussian quadrature formulas.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

A Gaussian quadrature rule is typically more accurate than a Newton-Cotes rule which requires the same number of function evaluations,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

if the integrand is smooth (i.e., if it has many derivatives.)

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Adaptive algorithms

If f does not have many derivatives at all points, or if the derivatives become large, then Gaussian quadrature is often insufficient. In this case, an algorithm similar to the following will perform better:

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

// This algorithm calculates the definite integral of a function

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

// from 0 to 1, adaptively, by choosing smaller steps near

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

// problematic points.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

// Set initial_h to the initial step size.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

x:=0

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

h:=initial_h

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

accumulator:=0

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

WHILE x

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

IF x+h>1.0 THEN

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

h=1.0-x

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

END IF

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

IF error from quadrature over for f is too large THEN

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Make h smaller

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

ELSE

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

accumulator:=accumulator + quadrature of f over

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

x:=x+h

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

IF error from quadrature over is very small THEN

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Make h larger

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

END IF

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

END IF

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

END WHILE

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

RETURN accumulator

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Some details of the algorithm require careful thought. For many cases, estimating the error from quadrature over an interval for a function f isn't obvious. One popular solution is to use two different rules of quadrature, and use their difference as an estimate of the error from quadrature. The other problem is deciding what "too large" or "very small" signify. A possible criterion for "too large" is that the quadrature error should not be larger than th where t, a real number, is the tolerance we wish to set for global error. Then again, if h is already tiny, it may not be worthwhile to make it even smaller even if the quadrature error is apparently large. This type of error analysis is usually called "a posteriori" since we compute the error after having computed the approximation.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Heuristics for adaptive quadrature are discussed by Forsythe et al. (Section 5.4).

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Extrapolation methods

The accuracy of a quadrature rule of the Newton-Cotes type is generally a function of the number of evaluation points.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

The result is usually more accurate as number of evaluation points increases,

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

or, equivalently, as the width of the step size between the points decreases.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

It is natural to ask what the result would be if the step size were allowed to approach zero.

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

This can be answered by extrapolating the result from two or more nonzero step sizes (see Richardson extrapolation).

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

The extrapolation function may be a polynomial or rational function.

Related Topics:
Polynomial - Rational function

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Extrapolation methods are described in more detail by Stoer and Bulirsch (Section 3.4).

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

Conservative (a priori) error estimation

Let f have a bounded first derivative over . The mean value theorem for f, where x < b, gives

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

: (x - a) f'(y_x) = f(x) - f(a),

~ ~ ~ ~ ~ ~ ~ ~ ~ ~

for some yx in depending on x. If we integrate in x from a to b on both sides and the take absolute values, we obtain

~ ~ ~ ~ ~ ~ ~ ~ ~ ~