Lecture 20: Quadrature

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

2.Gauss quadrature

Recall that the method of moments approach to numerical integration based upon sampling $𝒟=\{(x_i,y_i=f\left(x_i\right)),i=0,\dots ,n\}$,

imposes exact results for a finite number of members of a basis set $\{{\varphi }_0,\dots ,{\varphi }_n,\dots \}$

The trapezoid, Simpson formulas arise from the monomial basis set $\{1,t,t^2,\dots \}$, in which case

but any basis set can be chosen. Instead of prescribing the sampling points $x_i$ a priori, which typically leads to an error $e=𝒪\left({\varphi }_{n+1}\left(t\right)\right)$, the sampling points can be chosen to minimize the error $e$. For the monomial basis this leads to a system of $2(n+1)$ equations

for the unknown $n+1$ quadrature weights $w_i$ and the $n+1$ sampling points $x_i$. The system is nonlinear, but can be solved in an insightful manner exploiting the properties of orthogonal polynomials known as Gauss quadrature.

The basic idea is to consider a Hilbert function space with the scalar product

and orthonormal basis set $\{{\varphi }_0\left(t\right),{\varphi }_1\left(t\right),{\varphi }_2\left(t\right),,\dots \}$,

Assume that ${\varphi }_k\left(t\right)$ are polynomials of degree $k$. A polynomial $p_{2n+1}$ of degree $2n+1$ can be factored as

where $q_n\left(t\right)$ is the quotient polynomial of degree $n$, and $r_n$ is the remainder polynomial of degree $n$. The weighted integral of $p_{2n+1}$ is therefore

Since $\{{\varphi }_0,\dots ,{\varphi }_{n+1}\}$ is an orthonormal set, $(q_n,{\varphi }_{n+1})=0$, and the integral becomes

The integral of the $n^{\mathrm{th}}$ remainder polynomial can be exactly evaluated through an $n+1$ point quadrature

that however evaluates $r\left(t\right)$ rather than the original integrand $p_{2n+1}\left(t\right)$. However, evaluation of the factorization (2) at the roots $x_i$ of ${\varphi }_{n+1}$, ${\varphi }_{n+1}\left(x_i\right)=0$, $i=0,1,\dots ,n$, gives

stating that the values of the remainder at these nodes are the same as those of the $p_{2n+1}$ polynomial. This implies that

is an exact quadrature of order $2n+1$, $e=𝒪\left(t^{2n+1}\right)$. The weights $w_i$ can be determined through any of the previously outlined methods, e.g., method of moments

which is now a linear system that can be readily solved. Alternatively, the weights are also directly given as integrals of the Lagrange polynomials based upon the nodes that are roots of ${\varphi }_{n+1}$