Editing Numerical integration (section)

==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.