1.Numerical integration

Numerical integration proceeds along the same lines as numerical differentiation

with a different operator

with $ℒ:C\left([a,b]\backslash Z\right)\to ℝ$ with $Z$ a set of measure zero, e.g., $Z=\varnothing$ for a function continuous at all points in $[a,b]$. It is often useful to explicitly identify a weight function $\omega \left(t\right)$ that can attribute higher significance to subsets of $[a,b]$. In the integration case, the approximation basis choice can be combined with decomposition of the domain into $m$ subintervals $[a_k,b_k)$ such that ${\bigcup }_{k=1}^m[a_k,b_k)=[a,b)$ through

with

Monomial basis.

Integration of $f$ for $\omega \left(t\right)=1$ is approximated as

As for numerical differentation, the computational effort in the above formulation can be reduced through alternative basis choices.

Lagrange basis.

Trapezoid formula

For $d=1$, a linear approximant is defined over each subinterval $[x_{k-1},x_k]$, stated in Lagrange form as

${\displaystyle p_k\left(t\right)={\ell }_{k-1}\left(t\right)y_{k-1}+{\ell }_k\left(t\right)y_k=\frac{t-x_k}{x_{k-1}-x_k}{y}_{k-1}+\frac{t-x_{k-1}}{x_k-x_{k-1}}{y}_k.}$

The resulting integral approximation

${\displaystyle ℒ\left(f\right)=\underset{x_{k-1}}{\overset{x_k}{\int }}f\left(t\right)dt\cong ℐ\left(p\right)=\underset{x_{k-1}}{\overset{x_k}{\int }}{p}_k\left(t\right)dt=(x_k-x_{k-1})\frac{y_{k-1}+y_k}{2}.}$

The approximation error results from the known polynomial interpolation error formula

${\displaystyle f\left(t\right)-p\left(t\right)=\frac{f^{\text{'}\text{'}}\left({\xi }_t\right)}{2}(t-x_k)(t-x_{k-1}),}$

that gives

${\displaystyle e_k\left(t\right)=|\underset{x_{k-1}}{\overset{x_k}{\int }}[f\left(t\right)-p\left(t\right)]dt|⩽{||f^{\text{'}\text{'}}||}_{\infty }\frac{{(x_k-x_{k-1})}^3}{12}.}$

Using an equidistant partition $x_k=a+kh$, $h=(b-a)/n$, over the entire $[a,b]$ interval gives the composite trapezoid formula

${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt=\sum _{k=1}^n\underset{x_{k-1}}{\overset{x_k}{\int }}f\left(t\right)dt\cong \sum _{k=1}^n\underset{x_{k-1}}{\overset{x_k}{\int }}{p}_k\left(t\right)dt=h(\frac{y_0}{2}+\sum _{k=1}^n{y}_k+\frac{y_n}{2}).}$

The overall approximation error is bounded by the error over each subinterval

${\displaystyle e⩽\sum _{k=1}^n{e}_k⩽{||f^{\text{'}\text{'}}||}_{\infty }n\frac{h^3}{12}=(b-a)\frac{h^2{||f^{\text{'}\text{'}}||}_{\infty }}{12},}$

and the trapezoid formula exhibits quadratic convergence, $e=𝒪\left(h^2\right)$.

Simpson formula

A more accurate, quadratic approximant is obtained for $d=2$ using sample points $x_{k-2}$, $x_{k-1},x_k$

${\displaystyle p_k\left(t\right)={\ell }_{k-2}\left(t\right)y_{k-2}+{\ell }_{k-1}\left(t\right)y_{k-1}+{\ell }_k\left(t\right)y_k,}$ ${\displaystyle p_k\left(t\right)=\frac{(t-x_{k-1})(t-x_k)}{(x_{k-2}-x_{k-1})(x_{k-2}-x_k)}{y}_{k-2}+\frac{(t-x_{k-2})(t-x_k)}{(x_{k-1}-x_{k-2})(x_{k-1}-x_k)}{y}_{k-1}+\frac{(t-x_{k-2})(t-x_{k-1})}{(x_k-x_{k-2})(x_k-x_{k-1})}{y}_k.}$

Assuming $n=2m$, $x_k=a+kh$, $h=(b-a)/n$, over a subinterval the Simpson approximation is

${\displaystyle \underset{x_{k-2}}{\overset{x_k}{\int }}f\left(t\right)dt\cong \underset{x_{k-2}}{\overset{x_k}{\int }}{p}_k\left(t\right)dt=\frac{h}{3}(y_{k-2}+4y_{k-1}+y_k).}$

Integration of the interpolation error

${\displaystyle f\left(t\right)-p_k\left(t\right)=\frac{f^{\left(3\right)}\left({\xi }_t\right)}{3!}(t-x_{k-2})(t-x_{k-1})(t-x_k),}$

leads to calculation of

${\displaystyle \underset{x_{k-2}}{\overset{x_k}{\int }}(t-x_{k-2})(t-x_{k-1})(t-x_k)dt=\underset{-h}{\overset{h}{\int }}(\tau -h)\tau (\tau +h)d\tau =0.}$

(1)

This null result is a feature of even-degree interpolants $d=2r$. Note that the interpolation error formula contains an evaluation of $f^{\left(3\right)}\left({\xi }_t\right)$ at some point ${\xi }_t\in [x_{k-2},x_k]$ that depends on $t$, so the integral

${\displaystyle \underset{x_{k-2}}{\overset{x_k}{\int }}\frac{f^{\left(3\right)}\left({\xi }_t\right)}{3!}(t-x_{k-2})(t-x_{k-1})(t-x_k)dt}$

is not necessarily equal to zero. To obtain an error estimate, rewrite the interpolating polynomial in Newton form

${\displaystyle p_k\left(t\right)=y_{k-1}+[y_{k-1},y_{k-2}](t-x_{k-2})+[y_k,y_{k-1},y_{k-2}](t-x_{k-2})(t-x_{k-1}),}$

The next higher degree interpolant would be

${\displaystyle q_k\left(t\right)=p_k\left(t\right)+c_k(t-x_{k-2})(t-x_{k-1})(t-x_k),}$

and (1) implies that the integral of $q_k$ is equal to that of $p_k$

${\displaystyle \underset{x_{k-2}}{\overset{x_k}{\int }}{q}_k\left(t\right)dt=\underset{x_{k-2}}{\overset{x_k}{\int }}{p}_k\left(t\right)dt,}$

hence the Simpson formula based on a quadratic interpolation is as accurate as that based on a cubic interpolation. The error can now be estimated using

${\displaystyle e_k=|\underset{x_{k-2}}{\overset{x_k}{\int }}[f\left(t\right)-q_k\left(t\right)]dt|⩽\frac{{||f^{\left(\mathrm{iv}\right)}||}_{\infty }}{4!}\underset{x_{k-2}}{\overset{x_k}{\int }}(t-x_{k-2})(t-x_{k-1})(t-x_k)(t-z_k)dt,}$

where $z_k$ is some additional interpolation point within $[x_{k-2},x_k]$. It is convenient to choose $z_k=x_{k-1}$, which corresponds to a Hermite interpolation condition at the midpoint. This is simply for the purpose of obtaining an error estimate, and does not affect the Simpson estimate of the integral. Carrying out the calculations gives

${\displaystyle e_k⩽\frac{{||f^{\left(4\right)}||}_{\infty }}{90}{h}^5.}$

When applied to the overall interval $[a,b]$, the Simpson formula is stated as

${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong \frac{b-a}{6}[f\left(a\right)+4f\left(\frac{a+b}{2}\right)+f\left(b\right)],}$

with error

${\displaystyle e⩽\frac{{||f^{\left(4\right)}||}_{\infty }}{2880}{(b-a)}^5.}$

The composite Simpson formula is

${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt=\sum _{k=1}^m\underset{x_{2k-2}}{\overset{x_{2k}}{\int }}f\left(t\right)dt\cong \sum _{k=1}^m\underset{x_{2k-2}}{\overset{x_{2k}}{\int }}{p}_k\left(t\right)dt=\frac{h}{3}(y_0+4\sum _{k=1}^m{y}_{2k-1}+2\sum _{k=1}^{m-1}{y}_{2k}+y_n),}$

with an overall error bound

${\displaystyle e⩽\sum _{k=1}^m\underset{x_{2k-2}}{\overset{x_{2k}}{\int }}|f\left(t\right)-p_k\left(t\right)|dt⩽\frac{{||f^{\left(4\right)}||}_{\infty }}{90}{h}^5\frac{n}{2}=\frac{(b-a)}{180}{||f^{\left(4\right)}||}_{\infty }{h}^4,}$

that shows fourth-order convergence $e=𝒪\left(h^4\right)$.

Moment method.

using data $𝒟=\{(x_i,y_i=f\left(x_i\right)),i=0,\dots ,n\}$. Once the sampling points are chosen, there remain $n+1$ weights $w_i$ to be determined. One approach is to enforce exact quadrature for the first $n+1$ monomials

A Vandermonde system results whose solution furnishes the appropriate weights. Since the Vandermonde matrix is ill-conditioned, a solution is sought in exact arithmetic, using rational numbers as opposed to floating point approximations.

The principal utility of the moment method is construction of quadrature formulas for singular integrands. For example, in the integral

the $\ln t$ is an integrable singularity, and accurate quadrature rules can be constructed by the method of moments for

with $\omega \left(t\right)=-\ln t$.