MATH661

1.Gauss quadrature

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

\int_{a}^{b} ω (t) f (t) d t = \sum_{i = 0}^{n} w_{i} y_{i} + e ≅ \sum_{i = 0}^{n} w_{i} y_{i},

imposes exact results for a finite number of members of a basis set ${ϕ_{0}, \dots, ϕ_{n}, \dots}$

\int_{a}^{b} ω (t) ϕ_{k} (t) d t = \sum_{i = 0}^{n} w_{i} ϕ_{k} (x_{i}), k = 0, 1, \dots, n .

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

\int_{a}^{b} ω (t) t^{k} d t = \sum_{i = 0}^{n} w_{i} x_{i}^{k}, k = 0, 1, \dots, n,

but any basis set can be chosen. Instead of prescribing the sampling points $x_{i}$ a priori, which typically leads to an error $e = 𝒪 (ϕ_{n + 1} (t))$ , 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

\int_{a}^{b} ω (t) t^{k} d t = \sum_{i = 0}^{n} w_{i} x_{i}^{k}, k = 0, 1, \dots, 2 n + 1,

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

(f, g) = \int_{a}^{b} ω (t) f (t) g (t) d t,

and orthonormal basis set ${ϕ_{0} (t), ϕ_{1} (t), ϕ_{2} (t),, \dots}$ ,

(ϕ_{j}, ϕ_{k}) = \int_{a}^{b} ω (t) ϕ_{j} (t) ϕ_{k} (t) d t = δ_{j k} .

Assume that $ϕ_{k} (t)$ are polynomials of degree $k$ . A polynomial $p_{2 n + 1}$ of degree $2 n + 1$ can be factored as

p_{2 n + 1} (t) = q_{n} (t) ϕ_{n + 1} (t) + r_{n} (t),

where $q_{n} (t)$ is the quotient polynomial of degree $n$ , and $r_{n}$ is the remainder polynomial of degree $n$ . The weighted integral of $p_{2 n + 1}$ is therefore

\int_{a}^{b} ω (t) p_{2 n + 1} (t) d t = \int_{a}^{b} ω (t) [q_{n} (t) ϕ_{n + 1} (t) + r_{n} (t)] d t = (q_{n}, ϕ_{n + 1}) + \int_{a}^{b} ω (t) r_{n} (t) d t .

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

\int_{a}^{b} ω (t) p_{2 n + 1} (t) d t = \int_{a}^{b} ω (t) r_{n} (t) d t .

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

\int_{a}^{b} ω (t) r_{n} (t) d t = \sum_{i = 0}^{n} w_{i} r (x_{i}),

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

p_{2 n + 1} (x_{i}) = q_{n} (x_{i}) ϕ_{n + 1} (x_{i}) + r_{n} (x_{i}) = r_{n} (x_{i}),

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

\int_{a}^{b} ω (t) p_{2 n + 1} (t) d t = \sum_{i = 0}^{n} w_{i} p_{2 n + 1} (x_{i}),

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

\int_{a}^{b} ω (t) t^{k} d t = \sum_{i = 0}^{n} w_{i} x_{i}^{k}, k = 0, \dots, n,

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 $ϕ_{n + 1}$

w_{i} = \int_{a}^{b} ω (t) ℓ_{i} (t) d t .