Adaptive quadrature

Last updated

Adaptive quadrature is a numerical integration method in which the integral of a function is approximated using static quadrature rules on adaptively refined subintervals of the region of integration. Generally, adaptive algorithms are just as efficient and effective as traditional algorithms for "well behaved" integrands, but are also effective for "badly behaved" integrands for which traditional algorithms may fail.

Contents

General scheme

Adaptive quadrature follows the general scheme

1. procedure integrate ( f, a, b, τ ) 2.      3.      4.     ifε > τthen 5.         m = (a + b) / 2 6.         Q = integrate(f, a, m, τ/2) + integrate(f, m, b, τ/2) 7.     endif 8.     return Q

An approximation to the integral of over the interval is computed (line 2), as well as an error estimate (line 3). If the estimated error is larger than the required tolerance (line 4), the interval is subdivided (line 5) and the quadrature is applied on both halves separately (line 6). Either the initial estimate or the sum of the recursively computed halves is returned (line 7).

The important components are the quadrature rule itself

the error estimator

and the logic for deciding which interval to subdivide, and when to terminate.

There are several variants of this scheme. The most common will be discussed later.

Basic rules

The quadrature rules generally have the form

where the nodes and weights are generally precomputed.

In the simplest case, Newton–Cotes formulas of even degree are used, where the nodes are evenly spaced in the interval:

When such rules are used, the points at which has been evaluated can be re-used upon recursion:

Newton-Cotes re-use.png

A similar strategy is used with Clenshaw–Curtis quadrature, where the nodes are chosen as

Or, when Fejér quadrature is used,

Other quadrature rules, such as Gaussian quadrature or Gauss-Kronrod quadrature, may also be used.

An algorithm may elect to use different quadrature methods on different subintervals, for example using a high-order method only where the integrand is smooth.

Error estimation

Some quadrature algorithms generate a sequence of results which should approach the correct value. Otherwise one can use a "null rule" which has the form of the above quadrature rule, but whose value would be zero for a simple integrand (for example, if the integrand were a polynomial of the appropriate degree).

See:

Subdivision logic

"Local" adaptive quadrature makes the acceptable error for a given interval proportional to the length of that interval. This criterion can be difficult to satisfy if the integrands are badly behaved at only a few points, for example with a few step discontinuities. Alternatively, one could require only that the sum of the errors on each of the subintervals be less than the user's requirement. This would be "global" adaptive quadrature. Global adaptive quadrature can be more efficient (using fewer evaluations of the integrand) but is generally more complex to program and may require more working space to record information on the current set of intervals.

See also

Notes

    Related Research Articles

    <span class="mw-page-title-main">Antiderivative</span> Concept in calculus

    In calculus, an antiderivative, inverse derivative, primitive function, primitive integral or indefinite integral of a function f is a differentiable function F whose derivative is equal to the original function f. This can be stated symbolically as F' = f. The process of solving for antiderivatives is called antidifferentiation, and its opposite operation is called differentiation, which is the process of finding a derivative. Antiderivatives are often denoted by capital Roman letters such as F and G.

    <span class="mw-page-title-main">Integral</span> Operation in mathematical calculus

    In mathematics, an integral is the continuous analog of a sum, which is used to calculate areas, volumes, and their generalizations. Integration, the process of computing an integral, is one of the two fundamental operations of calculus, the other being differentiation. Integration started as a method to solve problems in mathematics and physics, such as finding the area under a curve, or determining displacement from velocity. Today integration is used in a wide variety of scientific fields.

    <span class="mw-page-title-main">Riemann integral</span> Basic integral in elementary calculus

    In the branch of mathematics known as real analysis, the Riemann integral, created by Bernhard Riemann, was the first rigorous definition of the integral of a function on an interval. It was presented to the faculty at the University of Göttingen in 1854, but not published in a journal until 1868. For many functions and practical applications, the Riemann integral can be evaluated by the fundamental theorem of calculus or approximated by numerical integration, or simulated using Monte Carlo Integration.

    <span class="mw-page-title-main">Dirac delta function</span> Generalized function whose value is zero everywhere except at zero

    In mathematical physics, the Dirac delta distribution, also known as the unit impulse, is a generalized function or distribution over the real numbers, whose value is zero everywhere except at zero, and whose integral over the entire real line is equal to one.

    <span class="mw-page-title-main">Gaussian quadrature</span> Approximation of the definite integral of a function

    In numerical analysis, an n-point Gaussian quadrature rule, named after Carl Friedrich Gauss, is a quadrature rule constructed to yield an exact result for polynomials of degree 2n − 1 or less by a suitable choice of the nodes xi and weights wi for i = 1, ..., n.

    <span class="mw-page-title-main">Numerical integration</span> Methods of calculating definite integrals

    In analysis, numerical integration comprises a broad family of algorithms for calculating the numerical value of a definite integral. The term numerical quadrature is more or less a synonym for "numerical integration", especially as applied to one-dimensional integrals. Some authors refer to numerical integration over more than one dimension as cubature; others take "quadrature" to include higher-dimensional integration.

    <span class="mw-page-title-main">Simpson's rule</span> Method for numerical integration

    In numerical integration, Simpson's rules are several approximations for definite integrals, named after Thomas Simpson (1710–1761).

    <span class="mw-page-title-main">Newton–Cotes formulas</span>

    In numerical analysis, the Newton–Cotes formulas, also called the Newton–Cotes quadrature rules or simply Newton–Cotes rules, are a group of formulas for numerical integration based on evaluating the integrand at equally spaced points. They are named after Isaac Newton and Roger Cotes.

    In the calculus of variations and classical mechanics, the Euler–Lagrange equations are a system of second-order ordinary differential equations whose solutions are stationary points of the given action functional. The equations were discovered in the 1750s by Swiss mathematician Leonhard Euler and Italian mathematician Joseph-Louis Lagrange.

    In mathematics, the Riemann–Stieltjes integral is a generalization of the Riemann integral, named after Bernhard Riemann and Thomas Joannes Stieltjes. The definition of this integral was first published in 1894 by Stieltjes. It serves as an instructive and useful precursor of the Lebesgue integral, and an invaluable tool in unifying equivalent forms of statistical theorems that apply to discrete and continuous probability.

    In mathematics, the Cauchy principal value, named after Augustin Louis Cauchy, is a method for assigning values to certain improper integrals which would otherwise be undefined. In this method, a singularity on an integral interval is avoided by limiting the integral interval to the singularity.

    <span class="mw-page-title-main">Path integral formulation</span> Formulation of quantum mechanics

    The path integral formulation is a description in quantum mechanics that generalizes the action principle of classical mechanics. It replaces the classical notion of a single, unique classical trajectory for a system with a sum, or functional integral, over an infinity of quantum-mechanically possible trajectories to compute a quantum amplitude.

    <span class="mw-page-title-main">Quasi-Monte Carlo method</span> Numerical integration process

    In numerical analysis, the quasi-Monte Carlo method is a method for numerical integration and solving some other problems using low-discrepancy sequences. This is in contrast to the regular Monte Carlo method or Monte Carlo integration, which are based on sequences of pseudorandom numbers.

    <span class="mw-page-title-main">Monte Carlo integration</span> Numerical technique

    In mathematics, Monte Carlo integration is a technique for numerical integration using random numbers. It is a particular Monte Carlo method that numerically computes a definite integral. While other algorithms usually evaluate the integrand at a regular grid, Monte Carlo randomly chooses points at which the integrand is evaluated. This method is particularly useful for higher-dimensional integrals.

    <span class="mw-page-title-main">Euler–Bernoulli beam theory</span> Method for load calculation in construction

    Euler–Bernoulli beam theory is a simplification of the linear theory of elasticity which provides a means of calculating the load-carrying and deflection characteristics of beams. It covers the case corresponding to small deflections of a beam that is subjected to lateral loads only. By ignoring the effects of shear deformation and rotatory inertia, it is thus a special case of Timoshenko–Ehrenfest beam theory. It was first enunciated circa 1750, but was not applied on a large scale until the development of the Eiffel Tower and the Ferris wheel in the late 19th century. Following these successful demonstrations, it quickly became a cornerstone of engineering and an enabler of the Second Industrial Revolution.

    Adaptive Simpson's method, also called adaptive Simpson's rule, is a method of numerical integration proposed by G.F. Kuncir in 1962. It is probably the first recursive adaptive algorithm for numerical integration to appear in print, although more modern adaptive methods based on Gauss–Kronrod quadrature and Clenshaw–Curtis quadrature are now generally preferred. Adaptive Simpson's method uses an estimate of the error we get from calculating a definite integral using Simpson's rule. If the error exceeds a user-specified tolerance, the algorithm calls for subdividing the interval of integration in two and applying adaptive Simpson's method to each subinterval in a recursive manner. The technique is usually much more efficient than composite Simpson's rule since it uses fewer function evaluations in places where the function is well-approximated by a cubic function.

    Clenshaw–Curtis quadrature and Fejér quadrature are methods for numerical integration, or "quadrature", that are based on an expansion of the integrand in terms of Chebyshev polynomials. Equivalently, they employ a change of variables and use a discrete cosine transform (DCT) approximation for the cosine series. Besides having fast-converging accuracy comparable to Gaussian quadrature rules, Clenshaw–Curtis quadrature naturally leads to nested quadrature rules, which is important for both adaptive quadrature and multidimensional quadrature (cubature).

    The Gauss–Kronrod quadrature formula is an adaptive method for numerical integration. It is a variant of Gaussian quadrature, in which the evaluation points are chosen so that an accurate approximation can be computed by re-using the information produced by the computation of a less accurate approximation. It is an example of what is called a nested quadrature rule: for the same set of function evaluation points, it has two quadrature rules, one higher order and one lower order. The difference between these two approximations is used to estimate the calculational error of the integration.

    A Sommerfeld expansion is an approximation method developed by Arnold Sommerfeld for a certain class of integrals which are common in condensed matter and statistical physics. Physically, the integrals represent statistical averages using the Fermi–Dirac distribution.

    QUADPACK is a FORTRAN 77 library for numerical integration of one-dimensional functions. It was included in the SLATEC Common Mathematical Library and is therefore in the public domain. The individual subprograms are also available on netlib.

    References