Mathematics Source Library
C & ASM


Numerical Integration

Numerical Integration

Numerical integration is concerned with developing algorithms to approximate the integral of a function f(x). The most commonly used algorithms are Newton-Cotes formulas, Romberg's method, Gaussian quadrature, and to lesser extents Hermite's formulas and certain adaptive techniques.

Caveats

All of the numerical integration techniques listed above assume that the integrand is continuous on the interval of integration. The error estimates generally require that the integrands are differentiable of whatever order is required so that the formula for the error estimate makes sense. Below is a check list which one should verify before using any of the numerical algorithms.

Integrand Check List.

  • Are there any singularities of the integrand in the interval of integration?
    • If there are, then can they be removed?
    • A function which has jump discontinuities can be integrated by splitting the interval of integration into subintervals on which the function is continuous, and then numerically integrate over each subinterval and add the results.
    • Using integration by parts certain types of singular integrals can be expressed as the sum of a function and a non-singular integral.
    • In other cases, an approximation of the integral in a small neighborhood of the singularity can be obtained by hand and the interval of integration can be split into three intervals in which numerical integration can be used on two of them.
  • Does the function oscillate over the region of integration? If it does, then make sure that the step size is chosen to be smaller than the wave length of the function. The interval of integration can also be split into subintervals in which each subinterval is half a wave length and the algorithm is applied to each subinterval.

Newton-Cotes Methods

Newton-Cotes methods approximate the integral of a function by summing a weighted combination of the function evaluated at equally-spaced points, nodes. If the endpoints of the interval of integration are excluded from the sum, the method is called an open Newton-Cotes method and if the endpoints are included the method is called a closed Newton-Cotes method.

Romberg's Method

Romberg's method is used to estimate the integral of a function on a closed and bounded interval. Classically, the method consists of successively applying the composite trapezoidal rule each time halving the length of the subintervals and using a linear combination of the resulting sequence of estimates to estimate the integral by successively deleting the low order error terms in the Euler-Maclaurin summation formula. The process terminates when the change of the estimate is within a preassigned tolerance within a preassigned number of successive estimates.

Gaussian Quadrature

Instead of evaluating the integrand at equally spaced nodes as in Newton-Cotes methods, Gaussian quadrature methods make a judicious choice of nodes so as to maximize the precision of the numerical integration relative to the number of integrand evaluations. The common Gaussian quadrature methods are: Gauss-Legendre used to integrate a function f(x) over a closed and bounded interval [a,b], Gauss_Laguerre used to integrate a function of the form f(x) exp( -x ) over the positive x-axis { x : x > 0 }, Gauss-Hermite used to integrate a function of the form f(x) exp( -x² ) over the entire x-axis, { x : -∞ < x < ∞ }, and Gauss-Chebyshev used to integrate a function of the form f(x) / sqrt( 1-x² ) over the interval [-1,1].

Hermite's Quadrature Methods

Hermite's quadrature uses Hermite interpolation to approximate the integral. Hermite quadrature uses not only the values of the integrand evaluated at equally spaced nodes but also the first derivatives of the integrand evaluated at the nodes. Because the contributions due to the first derivatives are telescoping sums only the first derivatives at the endpoints of the integration interval contribute to the integral approximation. Also there exists quadrature formulas using higher order derivatives sometimes called Birkhoff quadrature. The Hermite quadrature using the first order derivatives of the integrand evaluated at the endpoints of the integration interval is also called the corrected trapezoidal rule.

C source code is available for the following Hermite quadrature method:

Adaptive Quadrature

Adaptive methods are used to estimate the integral of a function on a closed and bounded interval. The basic idea is to take smaller subintervals in regions where the function changes significantly and larger subintervals in regions where the function doesn't change significantly. What is significant depends on the quadrature technique involved, for Simpson's rule, if the function is nearly cubic (or quadratic, or linear, or constant) then the function doesn't change significantly from being cubic and a large interval can be used.