Numerical integration ($I(f)={\int }_a^{b}f(x)dx$)

Numerical integration approximates a definite integral using finitely many sampled function values. For an integrand $f$ on a finite interval $[a,b]$, $a<b$, a quadrature rule has the formThe quadrature nodes $x_j$ tell us where to evaluate $f$; the quadrature weights $w_j$ turn those values into an approximation to the accumulated signed area. The exact definite integral depends on the integrand throughout the interval, whereas the quadrature rule sees only these samples. Integrating a chosen interpolant exactly still generally approximates the original integrand's integral.^[ref]

The fundamental theorem of calculus evaluates $I(f)$ as $F(b)-F(a)$ when an antiderivative $F$ is available. Numerical integration is useful when no convenient antiderivative is known, or each evaluation of the integrand requires an expensive computation. The task then becomes choosing informative quadrature nodes and controlling the quadrature error with as few evaluations as possible. Fixed data also permit a quadrature rule, but refinement cannot supply new information between observations.^[ref]

For an exact comparison, take the polynomial $f(x)=x^2$ on $[0,1]$. Its antiderivative $x^3/3$ gives $I(f)=1/3$. On a partition of an interval into $n\ge 1$ equal panels, the composite midpoint quadrature rule uses $n$ midpoint samples, each with quadrature weight $1/n$; the composite trapezoidal rule uses $n+1$ endpoint samples, with outer quadrature weights $1/(2n)$ and interior quadrature weights $1/n$. The existing derivations under these two quadrature rules giveIn particular, $M_1=1/4$ and $T_1=1/2$ already bracket $1/3$; two panels give $M_2=5/16$ and $T_2=3/8$. The convex function $x^2$ explains the opposite signs. A midpoint rectangle can cross the curve within its panel, so its total signed area, rather than pointwise ordering, determines the underestimate.

Define the signed approximation error as $E_n=Q_n(f)-I(f)$. In this example the absolute errors are $1/(12n^2)$ and $1/(6n^2)$: doubling $n$ reduces each by exactly four. A plot with both axes on a logarithmic scale therefore has slope $-2$, expressing order two with respect to panel refinement. This is an exact calculation for this polynomial, including $n=1$, rather than an observed rate extrapolated from a few computed answers.

For a general integrand, the composite midpoint quadrature rule and composite trapezoidal rule give second-order quadrature error bounds when $f$ has two continuous derivatives on $[a,b]$ and a known bound on $|f^{\prime \prime }|$. These hypotheses describe behavior between samples. A formal error formula containing an unknown derivative is not a numerical accuracy certificate or a justified sample count. Finite samples do not certify quadrature accuracy makes the missing information explicit. Errors in the evaluations and roundoff error must also be considered; noise sensitivity of a quadrature rule treats the former.

More elaborate choices exploit other information: Gaussian quadrature obtains high polynomial exactness from carefully chosen quadrature nodes, and Simpson's rule uses a quadratic interpolant. Adaptive numerical integration compares approximations and concentrates samples where they disagree, often saving evaluations. Such a difference is an error estimate, not a certified bound without further assumptions: both approximations can miss the same feature.^[ref]

For any finite set of quadrature nodes $x_1,\dots ,x_m$ in $[a,b]$, $a<b$, consider the polynomialEvery sample is zero, so every quadrature rule at these quadrature nodes gives $Q(p)=0$. Yet $p$ is nonnegative and strictly positive away from its finitely many polynomial roots; its continuity makes the integral positive. Consequently the smooth functions $0$ and $g(x)=p(x)/A$ have identical samples but definite integrals $0$ and $1$. Multiplying $g$ by any positive constant makes the difference arbitrarily large. This is an information limit, even with exact arithmetic; merely assuming smoothness without quantitative bounds does not remove it. The related degree ceiling for quadrature exactness uses the same squared nodal polynomial to limit polynomial exactness.

For example, the endpoint quadrature nodes $0,1$ give $p(x)=x^2(1-x{)}^2$, $A=1/30$, and $g(x)=30x^2(1-x{)}^2$. Both the zero function and $g$ produce zero under the one-panel trapezoidal rule, although $I(g)=1$. More generally, any estimate based only on their shared samples has absolute error at least $1/2$ for one of these two integrands, by the triangle inequality.

The obstruction also applies to a deterministic adaptive procedure that terminates after finitely many evaluations on the zero function: form $p$ using every quadrature node it visited. On $g$, the same zero responses lead to the same choices and answer. Certified quadrature error control therefore needs additional information, such as a known bound on relevant derivatives, rather than agreement among sampled approximations alone.

Integrate the five Lagrange interpolation polynomials at the equally spaced nodes to obtain the displayed rule. It is exact for polynomials through degree four, and symmetry also makes it exact in degree five. Symmetric endpoint, inner-node and midpoint weights $a,b,c$ satisfy $2a+2b+c=2$, $2a+b/2=2/3$, $2a+b/8=2/5$, giving $a=7/45$, $b=32/45$, $c=12/45$.

A quadrature rule approximates an integral by a finite weighted sum of sampled function values:The quadrature nodes $x_i$ specify where to sample, and the quadrature weights $w_i$ specify how much each sample contributes. The density $\rho (x)$ specifies the target integral, whereas the discrete quadrature weights belong to the approximation. Ordinary numerical integration has $\rho =1$. Gaussian quadrature and the midpoint quadrature rule choose interior quadrature nodes; the trapezoidal rule includes endpoints.^[ref]

For ordinary numerical integration, exactness on the constant function $1$ requires ${\sum }_i{w}_i=b-a$; exactness on $x$ additionally requires ${\sum }_i{w}_i{x}_i=(b^2-a^2)/2$. These two checks make the quadrature rule exact for every affine function. For example, both the midpoint quadrature rule and trapezoidal rule give $6$ for ${\int }_0^2(2x+1)dx$. Testing progressively higher polynomial degrees defines the degree of precision of a quadrature rule. Exactness describes a class of integrands, rather than promising exact results for every function.

Composite midpoint quadrature rules and composite trapezoidal rules apply a simple quadrature rule on every panel of a partition of an interval. This buys accuracy by shrinking the panels while keeping the elementary formula fixed. For continuous functions sampled exactly, their quadrature errors tend to zero as the largest panel width tends to zero. The familiar second-order bounds require stronger assumptions about derivatives. The signed approximation error used here is $Q(f)-I_{\rho }(f)$; some sources define quadrature error with the opposite sign, so the convention must accompany any formula.

A quadrature node is a point at which the integrand is sampled by a quadrature rule. Under the affine map $x=(a+b)/2+(b-a)t/2$ from $[-1,1]$ to $[a,b]$, a reference quadrature node $t_i$ becomes $x_i=(a+b)/2+(b-a)t_i/2$; for ordinary numerical integration, the corresponding quadrature weight is multiplied by $(b-a)/2$.

A quadrature weight multiplies a sampled integrand value in a quadrature rule. For ordinary numerical integration, its units match those of the integration variable. Positive quadrature weights that sum to $b-a$ make $Q(f)/(b-a)$ a weighted average of the samples; signed quadrature weights can magnify sample perturbations, as quantified by noise sensitivity of a quadrature rule.

The degree of precision of a quadrature rule is the largest integer $d$ such that $Q(p)=I_{\rho }(p)$ for every polynomial $p$ of polynomial degree at most $d$. By linearity, checking the monomials $1,x,\dots ,x^d$ suffices. Exactness for one higher-degree polynomial alone does not raise this degree: every lower polynomial degree must also be covered.

If the observed samples are $y_i=f(x_i)+{\epsilon }_i$, the resulting quadrature rule changes by $\Delta Q={\sum }_i{w}_i{\epsilon }_i$. Bounded sample perturbations $|{\epsilon }_i|\le \delta$ therefore satisfyFor ordinary numerical integration with nonnegative quadrature weights exact on constant functions, the bound is $(b-a)\delta$, independent of panel refinement. A common bias ${\epsilon }_i=\delta$ attains it: the quadrature error due to the function approximation may vanish while this sample error persists. Interpolating a fixed data set onto a finer grid likewise does not create new observations or remove its uncertainty.

With new samples whose errors are independent random variables of zero expectation and common variance ${\sigma }^2$, a different model givesFor equal panels this is ${\sigma }^2(b-a{)}^2/n$ for the composite midpoint quadrature rule and ${\sigma }^2{h}^2(n-1/2)$ for the composite trapezoidal rule. These variance reductions require the stated independence and unbiasedness assumptions; a common bias or correlated sample errors need not shrink. The two displayed perturbation formulas follow directly from the triangle inequality and the variance additivity for independent random variables, respectively.

Integrate the Lagrange polynomial through values at all equally spaced nodes $0,\dots ,s$, including both endpoints. The resulting quadrature rule is exact through degree $s$. For even $s$, the node polynomial ${\prod }_{j=0}^s(t-j)$ is odd about $s/2$, so its integral vanishes and exactness extends through degree $s+1$. The weights need not remain positive as the number of nodes increases.

Simpson's rule approximates a definite integral by integrating the quadratic polynomial obtained by polynomial interpolation at two endpoints and their midpoint. For an interval $[a,b]$ with $a<b$, write $m=(a+b)/2$ and $h=(b-a)/2$. The three quadrature nodes are $a,m,b$, and the quadrature rule isThus its quadrature weights are $h/3,4h/3,h/3$. Their sum is $2h=b-a$, as required for exactness on constant functions. This is the three-node Newton-Cotes closed quadrature; its degree of precision of a quadrature rule is three, even though the interpolating polynomial has polynomial degree at most two.^[ref]

For an original worked example, take the cubic polynomial $f(x)=x^3-2x+4$ on the interval $[-1,2]$. Here $m=1/2$ and $h=3/2$, so the function evaluations giveSimpson's rule givesAn antiderivative supplies an independent exact definite integral:The quadratic polynomial through the same quadrature nodes is $p_2(x)=\frac{3}{2}{x}^2-\frac{1}{2}x+3$. It differs from $f$ between the quadrature nodes, but its definite integral agrees exactly. The picture shows why: the signed polynomial interpolation error has two lobes with opposite integrals.

Three function evaluations therefore capture more information about an integral than about the whole function. When an antiderivative is unavailable, composite Simpson's rule repeats this inexpensive local construction over a partition of an interval. For a sufficiently smooth function, halving the node spacing reduces its truncation error by approximately a factor of sixteen.

Polynomial interpolation makes Simpson's rule exact for every polynomial of polynomial degree at most two. To understand the extra polynomial degree, centre the coordinate at the midpoint: $y=x-m$, with endpoints $-h,h$. For any cubic polynomial $p$, let $p_2$ be its interpolating quadratic polynomial. The polynomial $p-p_2$ vanishes at all three quadrature nodes, sowhere $A$ is the coefficient of $x^3$ in $p$. As a function of $y$, this is an odd function, whose integral over the symmetric interval $[-h,h]$ is zero. Hence integrating $p_2$ also integrates $p$ exactly. This is the symmetry mechanism in the nodal-polynomial criterion for quadrature exactness.

In the worked example, the polynomial difference factors asIts definite integrals on $[-1,1/2]$ and $[1/2,2]$ are respectively $81/64$ and $-81/64$. Their cancellation concerns signed integrals, rather than pointwise accuracy of polynomial interpolation.

The degree of precision of a quadrature rule stops at three: on the same interval, the monomial $x^4$ gives $S(x^4)=69/8$, whereas ${\int }_{-1}^2{x}^{4}dx=33/5$. Their difference, exact integral minus approximation, is $-81/40$.

Composite Simpson's rule applies Simpson's rule to consecutive pairs of equal subintervals. Choose a positive even integer $n\ge 2$, put $h=(b-a)/n$, and take quadrature nodes $x_j=a+jh$ for $j=0,\dots ,n$. There are $n$ subintervals, $n/2$ two-subinterval panels, and $n+1$ quadrature nodes. ThenThe second sum is empty when $n=2$, recovering one panel. The quadrature weights, before multiplication by $h/3$, are $1,4,2,4,\dots ,2,4,1$. Every interior even-indexed quadrature node is an endpoint of two panels, so its two endpoint contributions add to two; every odd-indexed quadrature node is a midpoint with weight four. An odd $n$ cannot be covered by these pairs. The displayed formula requires equal spacing.^[ref]

For the exponential function on the interval $[0,1]$, the exact integral is $I=e-1$. With $n=2$, the formula gives $S_2\approx 1.718861151876593$, an overestimate of about $5.793234175\times 10^{-4}$. With $n=128$, it gives $S_{128}\approx 1.718281828494607$, an overestimate of about $3.556140064\times 10^{-11}$. The figure shows the intermediate refinements and compares every computed absolute error with the error bound for composite Simpson's rule.

For this exponential function, a geometric series also gives a compact independent formula. With $q=e^h$ and $h=1/n$, summing the $n/2$ panel contributions yieldsA Taylor series at $h=0$ givesThe Big O notation describes the remaining terms as $h$ tends to zero. Consequently $(S_n-I)/(S_{2n}-I)$ tends to sixteen: the computed ratios increase from about $15.65169$ to $15.99965$. Fourth-order truncation error is also obtained by Richardson extrapolation of the composite trapezoidal rule.^[ref]

Suppose $f$ has a fourth derivative that is a continuous function on the interval $[a,b]$, and $|f^{(4)}(x)|\le M$ there. With $I={\int }_a^{b}f(x)dx$, the signed quadrature error here is $I-S_n$. Then, for some point $\xi \in (a,b)$,In particular, a strictly positive fourth derivative makes $I-S_n<0$: composite Simpson's rule overestimates the integral.^[ref]

The bound's scale follows by adding the panel errors. A panel of width $2h$ has signed quadrature error $-h^5{f}^{(4)}({\xi }_j)/90$ and absolute error at most $M{h}^5/90$. There are $(b-a)/(2h)$ panels, so the triangle inequality gives $M(b-a)h^4/180$. The intermediate value theorem for the continuous fourth derivative replaces the average of the panel values $f^{(4)}({\xi }_j)$ by a single value $f^{(4)}(\xi )$, giving the signed formula.

The underlying single-panel calculation can be seen through the Peano kernel theorem. On the interval $[0,2]$, define the raw Peano kernel $K(t)=L[(x-t{)}_{+}^3]$, where $L$ is integral minus Simpson's rule and $(x-t{)}_{+}=\max (x-t,0)$ is the positive part. Then $K(t)=-t^3(4-3t)/12$ for $0\le t\le 1$ and $K(t)=-(2-t{)}^3(3t-2)/12$ for $1\le t\le 2$. It is strictly negative inside the interval and has integral $-1/15$. ThusThe uniform norm $\Vert f^{(4)}{\Vert }_{\infty }$ bounds the fourth derivative. The polynomial $f(x)=x^4/24$ attains equality. Scaling to endpoint spacing $h$ multiplies the coefficient by $h^5$, giving the sharp single-panel bound $h^5\Vert f^{(4)}{\Vert }_{\infty }/90$. The normalized Peano kernel is $K/3!$, so the external factorial must be omitted when using that normalization.

For the exponential function on $[0,1]$, $f^{(4)}(x)=e^x$, so $M=e$ and the bound is $e/(180n^4)$. Its leading actual absolute error is $(e-1)/(180n^4)$, smaller than the bound because a uniform upper bound on the fourth derivative does not use its variation across the interval. The theorem bounds truncation error with exact function evaluations; sample errors and roundoff error are additional contributions, as described by noise sensitivity of a quadrature rule.

On a compact interval, positive quadrature rule weights whose sum is the interval length define uniformly bounded functionals on continuous functions. If the rules eventually integrate every fixed polynomial exactly, the Weierstrass approximation theorem proves convergence for every continuous function.

An $(n+1)$-node quadrature rule already exact through degree $n$ is exact through degree $n+1+k$ exactly when its nodal polynomialis orthogonal to every polynomial of degree at most $k$.

No $(n+1)$-node quadrature rule for a positive interval weight can be exact through degree $2n+2$, because it evaluates the nonnegative nodal square $Q_{n+1}^2$ as zero although its integral is positive.

If an $(n+1)$-node quadrature rule for a positive interval weight is exact through degree $2n$, every weight is positive. Indeed, applying the rule to the square of the corresponding degree-$n$ Lagrange cardinal polynomial isolates that weight.

If each $(n+1)$-node quadrature rule is exact through degree $n$ and has positive weights, then it converges on every continuous function. The Weierstrass approximation theorem reduces the error to the uniform polynomial-approximation error, while exactness on constants controls the sum of the weights.

A quadrature rule is exact through degree $d$ precisely when its node weights reproduce the first $d+1$ moments of the integration measure.