In
analysis
Analysis ( : analyses) is the process of breaking a complex topic or substance into smaller parts in order to gain a better understanding of it. The technique has been applied in the study of mathematics and logic since before Aristotle (3 ...
, numerical integration comprises 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 the
numerical solution of differential equations. This article focuses on calculation of definite integrals.
The term numerical quadrature (often abbreviated to
''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.
The basic problem in numerical integration is to compute an approximate solution to a definite integral
:
to a given degree of accuracy. If is a smooth function integrated over a small number of dimensions, and the domain of integration is bounded, there are many methods for approximating the integral to the desired precision.
Reasons for numerical integration
There are several reasons for carrying out numerical integration, as opposed to analytical integration by finding the
antiderivative:
# The integrand ''f''(''x'') may be known only at certain points, such as obtained by
sampling. Some
embedded systems and other computer applications may need numerical integration for this reason.
# A formula for the integrand may be known, but it may be difficult or impossible to find an antiderivative that is an
elementary function. An example of such an integrand is ''f''(''x'') = exp(−''x''
''2''), the antiderivative of which (the
error function, times a constant) cannot be written in
elementary form
In mathematics, an elementary function is a function of a single variable (typically real or complex) that is defined as taking sums, products, roots and compositions of finitely many polynomial, rational, trigonometric, hyperbolic, and ex ...
.
# It may be possible to find an antiderivative symbolically, but it may be easier to compute a numerical approximation than to compute the antiderivative. That may be the case if the antiderivative is given as an infinite series or product, or if its evaluation requires a
special function that is not available.
History
The term "numerical integration" first appears in 1915 in the publication ''A Course in Interpolation and Numeric Integration for the Mathematical Laboratory'' by
David Gibb.
Quadrature is a historical mathematical term that means calculating area. Quadrature problems have served as one of the main sources of
mathematical analysis.
Mathematicians of Ancient Greece, according to the
Pythagorean doctrine, understood calculation of
area as the process of constructing geometrically a
square having the same area (''squaring''). That is why the process was named quadrature. For example, a
quadrature of the circle,
Lune of Hippocrates,
The Quadrature of the Parabola
''Quadrature of the Parabola'' ( el, Τετραγωνισμὸς παραβολῆς) is a treatise on geometry, written by Archimedes in the 3rd century BC and addressed to his Alexandrian acquaintance Dositheus. It contains 24 propositions rega ...
. This construction must be performed only by means of
compass and straightedge.
The ancient Babylonians used the
trapezoidal rule to integrate the motion of
Jupiter
Jupiter is the fifth planet from the Sun and the largest in the Solar System. It is a gas giant with a mass more than two and a half times that of all the other planets in the Solar System combined, but slightly less than one-thousand ...
along the
ecliptic.
For a quadrature of a rectangle with the sides ''a'' and ''b'' it is necessary to construct a square with the side
(the
Geometric mean of ''a'' and ''b''). For this purpose it is possible to use the following fact: if we draw the circle with the sum of ''a'' and ''b'' as the diameter, then the height BH (from a point of their connection to crossing with a circle) equals their geometric mean. The similar geometrical construction solves a problem of a quadrature for a parallelogram and a triangle.
Problems of quadrature for curvilinear figures are much more difficult. The
quadrature of the circle with compass and straightedge had been proved in the 19th century to be impossible. Nevertheless, for some figures (for example the
Lune of Hippocrates) a quadrature can be performed. The quadratures of a sphere surface and a
parabola segment done by
Archimedes became the highest achievement of the antique analysis.
* The area of the surface of a sphere is equal to quadruple the area of a
great circle of this sphere.
* The area of a segment of the
parabola cut from it by a straight line is 4/3 the area of the triangle inscribed in this segment.
For the proof of the results Archimedes used the
Method of exhaustion of
Eudoxus.
In medieval Europe the quadrature meant calculation of area by any method. More often the
Method of indivisibles was used; it was less rigorous, but more simple and powerful. With its help
Galileo Galilei
Galileo di Vincenzo Bonaiuti de' Galilei (15 February 1564 – 8 January 1642) was an Italian astronomer, physicist and engineer, sometimes described as a polymath. Commonly referred to as Galileo, his name was pronounced (, ). He ...
and
Gilles de Roberval found the area of a
cycloid arch,
Grégoire de Saint-Vincent investigated the area under a
hyperbola (''Opus Geometricum'', 1647), and
Alphonse Antonio de Sarasa
Alphonse Antonio de Sarasa was a Jesuit mathematician who contributed to the understanding of logarithms, particularly as areas under a hyperbola.
Alphonse de Sarasa was born in 1618, in Nieuwpoort in Flanders. In 1632 he was admitted as a no ...
, de Saint-Vincent's pupil and commentator, noted the relation of this area to
logarithms.
John Wallis algebrised this method: he wrote in his ''Arithmetica Infinitorum'' (1656) series that we now call the
definite integral, and he calculated their values.
Isaac Barrow and
James Gregory made further progress: quadratures for some
algebraic curves
In mathematics, an affine algebraic plane curve is the zero set of a polynomial in two variables. A projective algebraic plane curve is the zero set in a projective plane of a homogeneous polynomial in three variables. An affine algebraic pl ...
and
spirals.
Christiaan Huygens successfully performed a quadrature of some
Solids of revolution.
The quadrature of the hyperbola by Saint-Vincent and de Sarasa provided a new
function, the
natural logarithm, of critical importance.
With the invention of
integral calculus came a universal method for area calculation. In response, the term quadrature has become traditional, and instead the modern phrase "''computation of a univariate definite integral''" is more common.
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. The integrand is evaluated at a finite set of points called integration points and a weighted sum of these values is used to approximate the integral. The integration points and weights depend on the specific method used and the accuracy required from the approximation.
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 that 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.
A 'brute force' kind of numerical integration can be done, if the integrand is reasonably well-behaved (i.e.
piecewise continuous
Continuity or continuous may refer to:
Mathematics
* Continuity (mathematics), the opposing concept to discreteness; common examples include
** Continuous probability distribution or random variable in probability and statistics
** Continuous g ...
and of
bounded variation), 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
In the mathematics, mathematical field of numerical analysis, interpolation is a type of estimation, a method of constructing (finding) new data points based on the range of a discrete set of known data points.
In engineering and science, one ...
functions that are easy to integrate. Typically these interpolating functions are
polynomials. In practice, since polynomials of very high degree tend to oscillate wildly, only polynomials of low degree are used, typically linear and quadratic.
The simplest method of this type is to let the interpolating function be a constant function (a polynomial of degree zero) that passes through the point
. This is called the ''midpoint rule'' or ''
rectangle rule''
The interpolating function may be a straight line (an
affine function, i.e. a polynomial of degree 1)
passing through the points
and
.
This is called the ''
trapezoidal rule''
For either one of these rules, we can make a more accurate approximation by breaking up the interval
into some number
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
where the subintervals have the form
with
and
Here we used subintervals of the same length
but one could also use intervals of varying length
.
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.
Quadrature rules with equally spaced points have the very convenient property of ''nesting''. The corresponding rule with each interval subdivided includes all the current points, so those integrand values can be re-used.
If we allow the intervals between interpolation points to vary, we find another group of quadrature formulas, such as the
Gaussian quadrature formulas. A Gaussian quadrature rule is typically more accurate than a Newton–Cotes rule that uses the same number of function evaluations, if the integrand is
smooth (i.e., if it is sufficiently differentiable). Other quadrature methods with varying intervals include
Clenshaw–Curtis quadrature (also called Fejér quadrature) methods, which do nest.
Gaussian quadrature rules do not nest, but the related
Gauss–Kronrod quadrature formulas do.
Generalized midpoint rule formula
A generalized midpoint rule formula is given by
or
where
denotes
-th derivative. For example, substituting
and
in the generalized midpoint rule formula, we obtain an equation of the inverse tangent
where
is
imaginary unit and
Since at each odd
the numerator of the integrand becomes
, the generalized midpoint rule formula can be reorganized as
The following example of Mathematica code generates the plot showing difference between inverse tangent and its approximation truncated at
and
:
f heta_, x_:= theta/(1 + theta^2 * x^2);
aTan heta_, M_, nMax_:=
2*Sum Function[x,_Evaluate[D[f[theta,_x_.html" ;"title=",_Evaluate[D[f[theta,_x.html" ;"title="Function[x, Evaluate[D[f[theta, x">Function[x, Evaluate[D[f[theta, x ">,_Evaluate[D[f[theta,_x.html" ;"title="Function[x, Evaluate[D[f[theta, x">Function[x, Evaluate[D[f[theta, x ][(m - 1/2)/
M])/((2*n + 1)!*(2*M)^(2*n + 1)), , ];
Plot[, ,
PlotRange -> All]
For a function
defined over interval
, its integral is
Therefore, we can apply the generalized midpoint integration formula above by assuming that
.
Adaptive algorithms
If ''f''(''x'') 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:
def calculate_definite_integral_of_f(f, initial_step_size):
"""
This algorithm calculates the definite integral of a function
from 0 to 1, adaptively, by choosing smaller steps near
problematic points.
"""
x = 0.0
h = initial_step_size
accumulator = 0.0
while x < 1.0:
if x + h > 1.0:
h = 1.0 - x # At end of unit interval, adjust last step to end at 1.
if error_too_big_in_quadrature_of_f_over_range(f, , x + h:
h = make_h_smaller(h)
else:
accumulator += quadrature_of_f_over_range(f, , x + h
x += h
if error_too_small_in_quadrature_of_over_range(f, , x + h:
h = make_h_larger(h) # Avoid wasting time on tiny steps.
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''(''x'') 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 ''local'' criterion for "too large" is that the quadrature error should not be larger than ''t'' ⋅ ''h'' 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. A ''global'' criterion is that the sum of errors on all the intervals should be less than ''t''. 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 the 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, using
series acceleration In mathematics, series acceleration is one of a collection of sequence transformations for improving the rate of convergence of a series. Techniques for series acceleration are often applied in numerical analysis, where they are used to improve th ...
methods such as
Richardson extrapolation
In numerical analysis, Richardson extrapolation is a sequence acceleration method used to improve the rate of convergence of a sequence of estimates of some value A^\ast = \lim_ A(h). In essence, given the value of A(h) for several values of h, ...
. The extrapolation function may be a
polynomial or
rational function. Extrapolation methods are described in more detail by Stoer and Bulirsch (Section 3.4) and are implemented in many of the routines in the
QUADPACK
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 ...
library.
Conservative (a priori) error estimation
Let
have a bounded first derivative over
i.e.
The
mean value theorem for
where
for_some_
gives
for some
gives
for some
depending on
.
If we integrate in
from
to
on both sides and take the absolute values, we obtain
We can further approximate the integral on the right-hand side by bringing the absolute value into the integrand, and replacing the term in
by an upper bound
where the
supremum
In mathematics, the infimum (abbreviated inf; plural infima) of a subset S of a partially ordered set P is a greatest element in P that is less than or equal to each element of S, if such an element exists. Consequently, the term ''greatest ...
was used to approximate.
Hence, if we approximate the integral
by the
quadrature rule our error is no greater than the right hand side of . We can convert this into an error analysis for the
Riemann sum, giving an upper bound of
for the error term of that particular approximation. (Note that this is precisely the error we calculated for the example
.) Using more derivatives, and by tweaking the quadrature, we can do a similar error analysis using a
Taylor series (using a partial sum with remainder term) for ''f''. This error analysis gives a strict upper bound on the error, if the derivatives of ''f'' are available.
This integration method can be combined with
interval arithmetic
Interval arithmetic (also known as interval mathematics, interval analysis, or interval computation) is a mathematical technique used to put bounds on rounding errors and measurement errors in mathematical computation. Numerical methods using ...
to produce
computer proofs and ''verified'' calculations.
Integrals over infinite intervals
Several methods exist for approximate integration over unbounded intervals. The standard technique involves specially derived quadrature rules, such as
Gauss-Hermite quadrature for integrals on the whole real line and
Gauss-Laguerre quadrature for integrals on the positive reals. Monte Carlo methods can also be used, or a change of variables to a finite interval; e.g., for the whole line one could use
and for semi-infinite intervals one could use
as possible transformations.
Multidimensional integrals
The quadrature rules discussed so far are all designed to compute one-dimensional integrals. To compute integrals in multiple dimensions, one approach is to phrase the multiple integral as repeated one-dimensional integrals by applying
Fubini's theorem (the tensor product rule). This approach requires the function evaluations to
grow exponentially as the number of dimensions increases. Three methods are known to overcome this so-called ''
curse of dimensionality''.
A great many additional techniques for forming multidimensional cubature integration rules for a variety of weighting functions are given in the monograph by Stroud.
Integration on the
sphere
A sphere () is a geometrical object that is a three-dimensional analogue to a two-dimensional circle. A sphere is the set of points that are all at the same distance from a given point in three-dimensional space.. That given point is the c ...
has been reviewed by Hesse et al. (2015).
[Kerstin Hesse, Ian H. Sloan, and Robert S. Womersley: Numerical Integration on the Sphere. In W. Freeden et al. (eds.), Handbook of Geomathematics, Springer: Berlin 2015, ]
Monte Carlo
Monte Carlo methods and
quasi-Monte Carlo methods are easy to apply to multi-dimensional integrals. They may yield greater accuracy for the same number of function evaluations than repeated integrations using one-dimensional methods.
A large class of useful Monte Carlo methods are the so-called
Markov chain Monte Carlo algorithms, which include the
Metropolis–Hastings algorithm
In statistics and statistical physics, the Metropolis–Hastings algorithm is a Markov chain Monte Carlo (MCMC) method for obtaining a sequence of random samples from a probability distribution from which direct sampling is difficult. This se ...
and
Gibbs sampling
In statistics, Gibbs sampling or a Gibbs sampler is a Markov chain Monte Carlo (MCMC) algorithm for obtaining a sequence of observations which are approximated from a specified multivariate probability distribution, when direct sampling is diff ...
.
Sparse grids
Sparse grids were originally developed by Smolyak for the quadrature of high-dimensional functions. The method is always based on a one-dimensional quadrature rule, but performs a more sophisticated combination of univariate results. However, whereas the tensor product rule guarantees that the weights of all of the cubature points will be positive if the weights of the quadrature points were positive, Smolyak's rule does not guarantee that the weights will all be positive.
Bayesian Quadrature
Bayesian quadrature is a statistical approach to the numerical problem of computing integrals and falls under the field of
probabilistic numerics. It can provide a full handling of the uncertainty over the solution of the integral expressed as a
Gaussian process posterior variance.
Connection with differential equations
The problem of evaluating the integral
:
can be reduced to an
initial value problem for an
ordinary differential equation by applying the first part of the
fundamental theorem of calculus. By differentiating both sides of the above with respect to the argument ''x'', it is seen that the function ''F'' satisfies
:
Methods developed for ordinary differential equations, such as
Runge–Kutta methods
In numerical analysis, the Runge–Kutta methods ( ) are a family of implicit and explicit iterative methods, which include the Euler method, used in temporal discretization for the approximate solutions of simultaneous nonlinear equations. T ...
, can be applied to the restated problem and thus be used to evaluate the integral. For instance, the standard fourth-order Runge–Kutta method applied to the differential equation yields Simpson's rule from above.
The differential equation
has a special form: the right-hand side contains only the independent variable (here
) and not the dependent variable (here
). This simplifies the theory and algorithms considerably. The problem of evaluating integrals is thus best studied in its own right.
See also
*
Numerical methods for ordinary differential equations
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 can also ...
*
Truncation error (numerical integration)
*
Clenshaw–Curtis quadrature
*
Gauss-Kronrod quadrature
*
Riemann Sum or
Riemann Integral
*
Trapezoidal rule
*
Romberg's method
*
Tanh-sinh quadrature
*
Nonelementary Integral
References
*
Philip J. Davis and
Philip Rabinowitz, ''Methods of Numerical Integration''.
*
George E. Forsythe, Michael A. Malcolm, and
Cleve B. Moler
Cleve Barry Moler is an American mathematician and computer programmer specializing in numerical analysis. In the mid to late 1970s, he was one of the authors of LINPACK and EISPACK, Fortran libraries for numerical computing. He invented MATLA ...
, ''Computer Methods for Mathematical Computations''. Englewood Cliffs, NJ: Prentice-Hall, 1977. ''(See Chapter 5.)''
*
*
Josef Stoer and
Roland Bulirsch
Roland Zdeněk Bulirsch (10 November 1932 – 21 September 2022) was a German mathematician specialising in numerical analysis. He studied and taught at the Technical University of Munich, and taught internationally as visiting professor. He was ...
, ''Introduction to Numerical Analysis''. New York: Springer-Verlag, 1980. ''(See Chapter 3.)''
*
Boyer, C. B., ''A History of Mathematics'', 2nd ed. rev. by
Uta C. Merzbach, New York: Wiley, 1989 (1991 pbk ed. ).
*
Eves, Howard, ''An Introduction to the History of Mathematics'', Saunders, 1990, ,
External links
Integration: Background, Simulations, etc.at Holistic Numerical Methods Institute
from Wolfram Mathworld
Lobatto quadrature formulafrom Encyclopedia of Mathematics
Implementations of many quadrature and cubature formulaewithin the free
Tracker Component Library
Tracker may refer to:
Arts and entertainment Fictional characters
* Tracker (''G.I. Joe''), in the ''G.I. Joe'' universe
* Tracker (''PAW Patrol''), in the animated television series ''PAW Patrol''
* Tracker Cameron, in the television series ' ...
.
SageMath Online Integrator
{{Authority control
Numerical analysis
*
Articles with example Python (programming language) code