View source for Numerical integration

{{Short description|Methods of calculating definite integrals}}
[[Image:Integral as region under curve.svg|thumb|Numerical integration is used to calculate a numerical approximation for the value <math>S</math>, the area under the curve defined by <math>f(x)</math>.]]
{{Differential equations}}

In [[Numerical analysis|analysis]], '''numerical integration''' comprises a broad family of [[algorithm]]s for calculating the numerical value of a definite [[integral]].
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''';<ref>{{MathWorld | urlname=Cubature | title=Cubature }}</ref> others take "quadrature" to include higher-dimensional integration.

The basic problem in numerical integration is to compute an approximate solution to a definite integral

:<math>\int_a^b f(x) \, dx</math>

to a given degree of accuracy. If {{math|''f''(''x'')}} 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.

Numerical integration has roots in the geometrical problem of finding a square with the same area as a given plane figure (''[[quadrature (geometry)|quadrature]]'' or ''squaring''), as in the [[quadrature of the circle]].
The term is also sometimes used to describe the [[numerical ordinary differential equations|numerical solution of differential equations]].

== Motivation and need ==
There are several reasons for carrying out numerical integration, as opposed to analytical integration by finding the [[antiderivative]]:

# The integrand {{math|''f'' (''x'')}} may be known only at certain points, such as obtained by [[sampling (statistics)|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 {{math|1=''f'' (''x'') = [[Exponential function|exp]](−''x''{{sup|2}})}}, the antiderivative of which (the [[error function]], times a constant) cannot be written in [[elementary form]] {{xref|(see also: [[Nonelementary integral]])}}.
# 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 ==
{{main|Quadrature (geometry)}}
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 (mathematician)|David Gibb]].<ref>{{cite web|url=http://jeff560.tripod.com/q.html|title=Earliest Known Uses of Some of the Words of Mathematics (Q)|website=jeff560.tripod.com|access-date=31 March 2018}}</ref>

"Quadrature" is a historical mathematical term that means calculating area. Quadrature problems have served as one of the main sources of [[mathematical analysis]]. [[Greek mathematics|Mathematicians of Ancient Greece]], according to the [[Pythagoreanism|Pythagorean]] doctrine, understood calculation of [[area]] as the process of constructing geometrically a [[square (geometry)|square]] having the same area {{nowrap|(''squaring''){{hsp}}{{mdash}}{{hsp}}}}that is why the process was named "quadrature". Examples include [[quadrature of the circle]], [[Lune of Hippocrates]], and the treatise {{em|[[Quadrature of the Parabola]]}}. This construction must be performed only by means of [[Straightedge and compass construction|compass and straightedge]].

The ancient Babylonians used the [[trapezoidal rule]] to integrate the motion of [[Jupiter (planet)|Jupiter]] along the [[ecliptic]].<ref>{{cite journal|author1=Mathieu Ossendrijver|title=Ancient Babylonian astronomers calculated Jupiter's position from the area under a time-velocity graph|journal=Science|date=Jan 29, 2016|doi=10.1126/science.aad8085|volume=351|issue=6272|pages=482–484|pmid=26823423|bibcode=2016Sci...351..482O|s2cid=206644971 }}</ref>

[[File:Geometric mean.svg|thumb|left|220px|Antique method of finding the [[Geometric mean]] ]]
{{pb}}{{np}}
For a quadrature of a rectangle with the sides ''a'' and ''b'' it is necessary to construct a square with the side <math>x =\sqrt {ab}</math> (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.

[[File:Parabola and inscribed triangle.svg|thumb|200px|{{center|The area of a segment of a parabola}}]]
{{pb}}{{np}}
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 [[The Quadrature of the Parabola|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.
{{np}}{{pb}}
To prove the results, Archimedes used the [[method of exhaustion]] of [[Eudoxus of Cnidus|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]] 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]], de Saint-Vincent's pupil and commentator, noted the relation of this area to [[logarithm]]s.

[[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 (mathematician)|James Gregory]] made further progress: quadratures for some [[algebraic curves]] and [[spiral]]s. [[Christiaan Huygens]] successfully performed a quadrature of some [[Solid of revolution|Solids of revolution]].

The quadrature of the hyperbola by Saint-Vincent and de Sarasa provided a new [[Function (mathematics)|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==
A '''quadrature rule''' is an approximation of the [[integral|definite integral]] of a [[function (mathematics)|function]], usually stated as a [[weighted sum]] of function values at specified points within the domain of integration.

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 error. Also, each evaluation takes time, and the integrand may be arbitrarily complicated.

===Quadrature rules based on step functions===
A "brute force" kind of numerical integration can be done, if the integrand is reasonably well-behaved (i.e. [[piecewise]] [[continuous function|continuous]] and of [[bounded variation]]), by evaluating the integrand with very small increments.

[[Image:Integration rectangle.svg|right|thumb|300px|Illustration of the rectangle rule.]]
This simplest method approximates the function by a [[step function]] (a piecewise constant function, or a segmented polynomial of degree zero) that passes through the point <math display="inline"> \left( \frac{a+b}{2}, f \left( \frac{a+b}{2} \right)\right) </math>. This is called the ''midpoint rule'' or ''[[rectangle method|rectangle rule]]''
<math display="block">\int_a^b f(x)\, dx \approx (b-a) f\left(\frac{a+b}{2}\right).</math>

===Quadrature rules based on interpolating functions===
A large class of quadrature rules can be derived by constructing [[interpolation|interpolating]] functions that are easy to integrate. Typically these interpolating functions are [[polynomial]]s.  In practice, since polynomials of very high degree tend to [[Runge's phenomenon|oscillate wildly]], only polynomials of low degree are used, typically linear and quadratic.

[[Image:Integration trapezoid.svg|right|thumb|300px|Illustration of the trapezoidal rule.]]
The interpolating function may be a straight line (an [[affine function]], i.e. a polynomial of degree 1)
passing through the points <math> \left( a, f(a)\right) </math> and <math> \left( b, f(b)\right) </math>.
This is called the ''[[trapezoidal rule]]''
<math display="block">\int_a^b f(x)\, dx \approx (b-a) \left(\frac{f(a) + f(b)}{2}\right).</math>

[[Image:Integration simpson.svg|right|thumb|300px|Illustration of Simpson's rule.]]
For either one of these rules, we can make a more accurate approximation by breaking up the interval <math> [a,b] </math> into some number <math> n </math> 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
<math display="block">\int_a^b f(x)\, dx \approx 
\frac{b-a}{n} \left( {f(a) \over 2} + \sum_{k=1}^{n-1} \left( f \left( a + k \frac{b-a}{n} \right) \right) + {f(b) \over 2} \right),</math>

where the subintervals have the form <math> [a+k h,a+ (k+1)h] \subset [a,b], </math> with <math display="inline">h = \frac{b - a}{n}</math> and <math>k = 0,\ldots,n-1. </math> Here we used subintervals of the same length <math> h </math> but one could also use intervals of varying length <math> \left( h_k \right)_k </math>.

Interpolation with polynomials evaluated at equally spaced points in <math> [a,b] </math> 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 function|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 formula]]s do.

=== Adaptive algorithms ===
{{excerpt|Adaptive quadrature}}

=== Extrapolation methods ===
The accuracy of a quadrature rule of the [[Newton–Cotes formulas|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]] methods such as [[Richardson extrapolation]]. 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]] library.

=== Conservative (a priori) error estimation ===
Let <math>f</math> have a bounded first derivative over <math>[a,b],</math> i.e. <math>f \in C^1([a,b]).</math> The [[mean value theorem]] for <math> f,</math> where <math>x \in [a,b),</math> gives
<math display="block"> (x - a) f'(\xi_x) = f(x) - f(a), </math>
for some <math> \xi_x \in (a,x] </math> depending on <math> x </math>.

If we integrate in <math> x </math> from <math> a </math> to <math> b </math> on both sides and take the absolute values, we obtain
<math display="block">
    \left| \int_a^b f(x)\, dx - (b - a) f(a) \right|
  = \left| \int_a^b (x - a) f'(\xi_x)\, dx \right| .
</math>

We can further approximate the integral on the right-hand side by bringing the absolute value into the integrand, and replacing the term in <math> f' </math> by an upper bound
{{NumBlk|:|
<math>
  \left| \int_a^b f(x)\, dx - (b - a) f(a) \right| 
  \leq {(b - a)^2 \over 2} \sup_{a \leq x \leq b} \left| f'(x) \right| ,
</math>
|{{EquationRef|1}}
}}
where the [[supremum]] was used to approximate.

Hence, if we approximate the integral <math display="inline"> \int_a^b f(x) \, dx </math> by the [[#Methods for one-dimensional integrals|quadrature rule]] <math> (b - a) f(a) </math> our error is no greater than the right hand side of {{EquationNote|1}}. We can convert this into an error analysis for the [[Riemann sum#Definition|Riemann sum]], giving an upper bound of
<math display="block">\frac{n^{-1}}{2} \sup_{0 \leq x \leq 1} \left| f'(x) \right|</math>
for the error term of that particular approximation. (Note that this is precisely the error we calculated for the example <math>f(x) = x</math>.) 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]] to produce [[computer proof]]s 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.<ref>{{cite book |last=Leader |first=Jeffery J. | author-link=Jeffery J. Leader| title=Numerical Analysis and Scientific Computation |year=2004 |publisher=Addison Wesley |isbn= 978-0-201-73499-7}}</ref> 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
<math display="block">
 \int_{-\infty}^{\infty} f(x) \, dx
= \int_{-1}^{+1} f\left( \frac{t}{1-t^2} \right) \frac{1+t^2}{\left(1-t^2\right)^2} \, dt,
</math>
and for semi-infinite intervals one could use
<math display="block">\begin{align}
\int_a^{\infty} f(x) \, dx &= \int_0^1 f\left(a + \frac{t}{1-t}\right) \frac{dt}{(1-t)^2}, \\
\int_{-\infty}^a f(x) \, dx &= \int_0^1 f\left(a - \frac{1-t}{t}\right) \frac{dt}{t^2},
\end{align}</math>
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 [[exponential growth|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.<ref name="StroudBook">{{cite book |last1=Stroud |first1=A. H. |title=Approximate Calculation of Multiple Integrals |url=https://archive.org/details/approximatecalcu0000stro_b8j7 |url-access=registration |date=1971 |publisher=Prentice-Hall Inc. |location=Cliffs, NJ|isbn=9780130438935 }}</ref>
Integration on the [[sphere]] has been reviewed by Hesse et al. (2015).<ref>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, {{doi|10.1007/978-3-642-54551-1_40}}</ref>

=== Monte Carlo ===
{{main|Monte Carlo integration}}

[[Monte Carlo method]]s and [[quasi-Monte Carlo method]]s 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.{{Citation needed|date=November 2018}}

A large class of useful Monte Carlo methods are the so-called [[Markov chain Monte Carlo]] algorithms, which include the [[Metropolis–Hastings algorithm]] and [[Gibbs sampling]].

=== Sparse grids ===
[[Sparse grid]]s 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 definite integral
:<math>F(x) = \int_a^x f(u)\, du</math>
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
:<math> \frac{d F(x)}{d x} = f(x), \quad F(a) = 0. </math>

[[Numerical methods for ordinary differential equations]], such as [[Runge–Kutta methods]], 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 <math>F'(x) = f(x)</math> has a special form: the right-hand side contains only the independent variable (here <math>x</math>) and not the dependent variable (here <math>F</math>). This simplifies the theory and algorithms considerably. The problem of evaluating integrals is thus best studied in its own right.

Conversely, the term "quadrature" may also be used for the solution of differential equations: "[[Linear_differential_equation#Types_of_solution|solving by quadrature]]" or "[[Ordinary_differential_equation#Reduction_to_quadratures|reduction to quadrature]]" means expressing its solution in terms of [[antiderivative|integrals]].

==See 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==
{{Reflist}}
* [[Philip J. Davis]] and [[Philip Rabinowitz (mathematician)|Philip Rabinowitz]], ''Methods of Numerical Integration''.
* [[George E. Forsythe]], Michael A. Malcolm, and [[Cleve B. Moler]], ''Computer Methods for Mathematical Computations''. Englewood Cliffs, NJ: Prentice-Hall, 1977. ''(See Chapter 5.)''
* {{Citation |last1=Press |author-link=William H. Press |first1=W.H.|last2=Teukolsky|author2-link=Saul Teukolsky|first2=S.A.|last3=Vetterling|first3=W.T.|last4=Flannery|first4=B.P.|year=2007|title=Numerical Recipes: The Art of Scientific Computing|edition=3rd|publisher=Cambridge University Press| location=New York|isbn=978-0-521-88068-8|chapter=Chapter 4. Integration of Functions|chapter-url=http://apps.nrbook.com/empanel/index.html?pg=155}}
* [[Josef Stoer]] and [[Roland Bulirsch]], ''Introduction to Numerical Analysis''. New York: Springer-Verlag, 1980. ''(See Chapter 3.)''
* [[Carl Benjamin Boyer|Boyer, C. B.]], ''A History of Mathematics'', 2nd ed. rev. by [[Uta Merzbach|Uta C. Merzbach]], New York: Wiley, 1989 {{isbn|0-471-09763-2}} (1991 pbk ed. {{isbn|0-471-54397-7}}).
* [[Howard Eves|Eves, Howard]], ''An Introduction to the History of Mathematics'', Saunders, 1990, {{isbn|0-03-029558-0}},
* S.L.Sobolev and V.L.Vaskevich: ''The Theory of Cubature Formulas'', Kluwer Academic, ISBN 0-7923-4631-9 (1997).

==External links==
{{commons category}}
* [http://numericalmethods.eng.usf.edu/mws/gen/07int/index.html Integration: Background, Simulations, etc.] at Holistic Numerical Methods Institute
* [http://mathworld.wolfram.com/LobattoQuadrature.html Lobatto Quadrature] from Wolfram Mathworld
* [https://www.encyclopediaofmath.org/index.php/Lobatto_quadrature_formula Lobatto quadrature formula] from Encyclopedia of Mathematics
* [https://github.com/USNavalResearchLaboratory/TrackerComponentLibrary/tree/master/Mathematical%20Functions/Numerical%20Integration/Cubature%20Points Implementations of many quadrature and cubature formulae] within the free [[Tracker Component Library]].
* [https://analyticphysics.com/Coding%20Methods/SageMath%20Online%20Integrator.htm SageMath Online Integrator]

{{Numerical integration}}
{{Differential equations topics}}
{{Authority control}}

[[Category:Numerical analysis]]
[[Category:Numerical integration| ]]