Gausiaan Quadrature

• The abscissas
  • The abscissas are not equally spaced

• The coefficients of integrand includes weights

• $[a,b]$,
  • interval
• $\omega (x)$,
  • nonnegative weight function on the interval $[a,b]$,

Assumption

• (a) $\{x\in [a,b]\mid \omega (x)\ge 0\}$ is measurable
• (b) all momemnts exists and finite

$$\begin{array}{}\text{(1)}& \mathrm{\forall }k=0,1,\dots ,{\mu }_k:={\int }_a^b{x}^k\omega (x)dx<0.\end{array}$$

• (c) For polynomials $s(x)$ which are nonnegative on $[a,b]$,

$${\int }_a^b\omega (x)s(x)dx=0\Rightarrow s(x)\equiv 0.$$

■

$$\begin{array}{}\text{(2)}& I(f):={\int }_a^b\omega (x)f(x)dx,\end{array}$$ $$\begin{array}{}\text{(3)}& \tilde{I}(f):=\sum _{i=1}^n{w}_{i}f(x_i).\end{array}$$ $$\begin{array}{rcl}{\overline{\mathrm{\Pi }}}_j& :=& \{p\mid p(x)=x^j+a_1{x}^{j-1}+\cdots +a_j,a_k\in \mathbb{R}\}\\ {\mathrm{\Pi }}_j& :=& \{p\mid \mathrm{degree}(p)\le j\}.\end{array}$$ $$\begin{array}{rcl}(f,g)& :=& {\int }_a^b\omega (x)f(x)g(x)dx\\ \text{(4)}& (f,f)& =& {\int }_a^b\omega (x)f(x{)}^{2}dx.\end{array}$$

Let $h$ be function.

$$\begin{array}{rcl}(hf,g)& =& {\int }_a^b\omega (x)h(x)f(x)g(x)dx\\ & =& {\int }_a^b\omega (x)f(x)h(x)g(x)dx\\ & =& (f,hg).\end{array}$$

Theorem 3.6.3

There exists $p_j\in {\overline{\mathrm{\Pi }}}_j$ for $j=0,1,\dots ,$ such that

$$\begin{array}{}\text{(5)}& i\ne k,(p_i,p_k)=0.\end{array}$$

These polynomials are uniquely defined by the recursions

$$\begin{array}{rcl}\text{(6)}& p_0(x)& \equiv & 1\text{(7)}& p_{i+1}(x)& =& (x-{\delta }_{i+1})p_i(x)-{\gamma }_{i+1}^2{p}_{i-1}(x)i\ge 0,{\delta }_{i+1}& =& \frac{(x{p}_i,p_i)}{(p_i,p_i)}i\ge 0,\\ \text{(8)}& {\gamma }_{i+1}^2& =& \{\begin{array}{ll}1& i=0\\ \frac{(p_i,p_i)}{(p_{i-1},p_{i-1})}& i\ge 1\end{array}.\end{array}$$

proof

We prove by induction. If $i=0$, the statements are true. Let us assume that the statments holds true up to $j\le i$. We’ll show there exists a unique polynomials $p_{i+1}\in {\overline{\mathrm{\Pi }}}_{i+1}$ with

$$\begin{array}{}\text{(9)}& (p_{i+1},p_j)=0j\le i.\end{array}$$

Let $p_{i+1}\in {\overline{\mathrm{\Pi }}}_{i+1}$. There exist constants $\delta ,c_0,\dots ,c_{i-1}$ such that

$$p_{i+1}(x)=x{p}_i(x)-\delta p_i(x)+\sum _{k=0}^{i-1}{c}_k{p}_k(x).$$

We have

$$\begin{array}{}\text{(10)}& (p_{i+1},p_j)& =& (x{p}_i,p_j)-\delta (p_i,p_j)+\sum _{k=0}^{i-1}{c}_k(p_k,p_j).\end{array}$$

The constants must satisfy

$$\begin{array}{}\text{(11)}& (p_{i+1},p_i)& =& (x{p}_i,p_j)-\delta (p_i,p_j)=0\text{(12)}& j<i,(p_{i+1},p_j)& =& (x{p}_i,p_j)+c_j(p_j,p_j)=0.\end{array}$$

From $\text{(???)}$, $(p_k,p_k)\ne 0$. In fact, suppose that $(p_k,p_k)=0$.

$$\begin{array}{}\text{(13)}& j\le i,(p_j,p_j)& =& {\int }_a^b\omega (x)p_j^2(x)dx& =& 0.\end{array}$$

Hence $p_j\equiv 0$. This contradicts to the assumption of the induction.

The euqation $\text{(11)}$ can be solved uniquely

$$\begin{array}{rcl}\delta & =& \frac{(x{p}_i,p_i)}{(p_i,p_i)}.\end{array}$$

By the induction hypothesis, for $j<i$,

$$\begin{array}{rcl}& & p_j(x)=(x-{\delta }_j)p_{j-1}(x)-{\gamma }_j^2{p}_{j-2}(x)\\ \text{(14)}& & \iff & x{p}_{j-1}(x)=p_j(x)+{\gamma }_j^2{p}_{j-2}(x)+{\delta }_j{p}_{j-1}(x)\end{array}$$ $$\begin{array}{rcl}{c}_j& =& -\frac{(x{p}_j,p_i)}{(p_j,p_j)}\\ & =& -\frac{(p_{j+1},p_i)+{\gamma }_{j+1}^2(p_{j-1},p_i)+{\delta }_{j+1}(p_j,p_i)}{(p_j,p_j)}\\ & =& -\frac{(p_{j+1},p_i)}{(p_j,p_j)}\\ & =& \{\begin{array}{ll}0& (j<i-1)\\ -\frac{(p_j,p_i)}{(p_j,p_j)}& (j=i-1)\end{array}.\end{array}$$

$◻$

Corollary 3.6.9

• $\{p_0,\dots ,p_n\}$,
  • polynomials in theorem 3.6.3

$$(p,p_n)=0p\in {\mathrm{\Pi }}_{n-1}.$$

proof

Since $p_j\in {\mathrm{\Pi }}_j$ for all $j$, $p$ is representable as a linear combination of the polynomials.

$◻$

Theorem 3.6.10

• $\{x_1,\dots x_n\}$,
  • the roots of $p_n$,

The roots are real and simple. $x_i\in (a,b)$.

proof

Assume that $a<x_1,\dots ,x_l<b$ are distinct roots of $p_n$ and $l<n$. Let

$$q(x):=\prod _{j=1}^l(x-x_j)\in {\overline{\mathrm{\Pi }}}_l.$$

Apparently, $q{p}_n\not\equiv 0$. That implies

$$\begin{array}{rcl}0& =& (p_n,q)(\because l<n\text{ and Theorem 3.6.3})\end{array}$$

On the other hands,

$$\begin{array}{rcl}(p_n,q)& =& {\int }_a^b\omega (x)p_n(x)q(x)dx\\ \text{(15)}& & \ne & 0(\because \text{Assumption (c)}).\end{array}$$

$◻$

Theorem 3.6.11

• $t_1,\dots ,t_n\in \mathbb{R}$,
  • $t_i\ne t_j$ for $i\ne j$,

$$A:=\left(\begin{array}{ccc}{p}_0(t_1)& \cdots & p_0(t_n)\\ ⋮& & ⋮\\ p_{n-1}(t_1)& \cdots & p_{n-1}(t_n)\end{array}\right).$$

$A$ is nonsingular.

proof

Assume that $A$ is singular. There exists a vector $c^{\mathrm{T}}=(c_0,\dots ,c_{n-1})\ne 0$ such that

$$c^{\mathrm{T}}A=0.$$

Let

$$q(x):=\sum _{i=0}^{n-1}{c}_i{p}_i(x).$$

Since

$$\begin{array}{rcl}{c}^{\mathrm{T}}A& =& (\sum _{i=0}^{n-1}{c}_i{p}_i(t_1),\cdots ,\sum _{i=0}^{n-1}{c}_i{p}_i(t_1))\\ & =& (q(t_1)\cdots ,q(t_n))\\ & =& (0,\dots ,0),\end{array}$$

$t_1,\dots ,t_n$ are distinct roots of $q$. Since $\mathrm{degree}(p_i)<n$ for all $i=0,\dots ,n-1$, $\mathrm{degree}(q)<n$. Thus, $q\equiv 0$. However, $\{p_i\}$ is linearly independet, so $c_i=0$ for all $i$.

$◻$

Definition Haar condition

• $\{p_i\}$,
  • polynomials

$\{p_i\}$ is said to satisfy Haar condition if for all distincts $t_1,\dots ,t_n$,

$$\left(\begin{array}{ccc}{p}_0(t_1)& \cdots & p_0(t_n)\\ ⋮& & ⋮\\ p_{n-1}(t_1)& \cdots & p_{n-1}(t_n)\end{array}\right)$$

is nonsingular. Such $\{p_i\}$ is said to be Chebychev system.

■

Theorem 3.6.12

• $x_1,\dots ,x_n$,
  • roots of $n$-th orthonormal polynomials $p_n$,
• $w_1,\dots ,w_n$,
  • solution of nonsingular system of equations

(a)

$$\begin{array}{}\text{(16)}& \left(\begin{array}{ccc}{p}_0(x_1)& & p_0(x_n)\\ & & \\ p_{n-1}(x_1)& & p_{n-1}(x_n)\end{array}\right)w& =& \left(\begin{array}{c}(p_0,p_0)\\ 0\\ ⋮\\ 0\end{array}\right)\end{array}$$

Then $w_i>0$ for all $i=1,\dots ,n$ and for all $p\in {\mathrm{\Pi }}_{2n-1}$

$$\begin{array}{}\text{(17)}& {\int }_a^b\omega (x)p(x)dx& =& \sum _{i=1}^n{w}_{i}p(x_i).\end{array}$$

(b)

Conversely, if $w_i,x_i$ satisfy $\text{(17)}$, then $x_i$ are the roots of $p_n$ and the weights $w_i$ satisfy $\text{(16)}$.

(c)

There is no $x_i,w_i,i=1,\dots ,n$ such that $\text{(17)}$ holds for all $p\in {\mathrm{\Pi }}_{2n}$.

proof

(a)

By Theorem 3.6.10, the roots $x_i,i=1,\dots ,n$ of $p_n$ are real and mutually distinct numbers in $(a,b)$. The matrix

$$\begin{array}{}\text{(18)}& A& :=& \left(\begin{array}{ccc}{p}_0(x_1)& \cdots & p_0(x_n)\\ ⋮& & ⋮\\ p_{n-1}(x_1)& \cdots & p_{n-1}(x_n)\end{array}\right)\end{array}$$

is nonsingular by Theorem 3.6.11. Hence $\text{(16)}$ has a unique solution. Let $p\in \prod _{2n-1}$. There exist $q,r\in {\mathrm{\Pi }}_{n-1}$ such that

$$p(x)=p_n(x)q(x)+r(x).$$

$p,q$ can be written in

$$\begin{array}{rcl}q(x)& =& \sum _{k=0}^{n-1}{\alpha }_k{p}_k(x)\\ \text{(19)}& r(x)& =& \sum _{k=0}^{n-1}{\alpha }_k{p}_k(x).\end{array}$$

Since $p_0\equiv 1$,

$$\begin{array}{rcl}{\int }_a^b\omega (x)p(x)dx& =& {\int }_a^b\omega (x)(p_n(x)q(x)+r(x)p_0(x))dx\\ & =& (p_n,q)+(r,p_0)\\ & =& (p_n,\sum _{k=0}^{n-1}{\alpha }_k{p}_k)+(\sum _{k=0}^{n-1}{\beta }_k{p}_k,p_0)\\ & =& {\beta }_0(p_0,p_0).\end{array}$$

On the other hand,

$$\begin{array}{rcl}\sum _{k=1}^n{w}_{i}p(x_i)& =& \sum _{i=1}^n{w}_i(p_n(x_i)q(x_i)+r(x_i))\\ & =& \sum _{k=0}^{n-1}{\beta }_k\sum _{i=1}^n{w}_i{p}_k(x_i)\\ & =& {\beta }_0\sum _{i=1}^n{w}_i{p}_0(x_i)\\ & =& {\beta }_0(p_0,p_0)\end{array}$$

We claim that if $w_i,x_i$, $i=1,\dots ,n$ are such that $\text{(17)}$ holds for all $p\in {\mathrm{\Pi }}_{2n-1}$, then $w_i>0$ for $i=1,\dots ,n$. Let

$${\overline{p}}_j(x):=\prod _{h=1,h\ne j}^n(x-x_h{)}^2\in {\mathrm{\Pi }}_{2n-2},(j=1,\dots ,n).$$

Since ${\overline{p}}_j\not\equiv 0$, we have by assumption (c),

$$\begin{array}{rcl}0& <& {\int }_a^{b}w(x){\overline{p}}_j(x)dx\\ & =& \sum _{i=1}^n{w}_i{\overline{p}}_j(x_i)\\ & =& w_j\prod _{h=1,h\ne j}^n(x_j-x_j{)}^2.\end{array}$$

Thus, $w_j>0$.

(c)

Assume that $w_i,x_i,i=1,\dots ,n$ satisfy $\text{(17)}$ for all $p\in {\mathrm{\Pi }}_{2n}$. Let

$$\overline{p}(x):=\prod _{j=1}^n(x-x_j{)}^2\in {\mathrm{\Pi }}_{2n}.$$ $$\begin{array}{rcl}0& <& {\int }_a^b\omega (x)\overline{p}(x)dx(\because \text{Assumption (c)})\\ & =& \sum _{i=1}^n{w}_i\overline{p}(x_i)\\ & =& 0.\end{array}$$

(b)

Suppose that $w_i,x_i$ are such that $\text{(17)}$ for all $p\in {\mathrm{\Pi }}_{2n-1}$. $x_i$ must be $x_i\ne x_{j}i\ne j$. Indeed, let us assume $x_{n-1}=x_n$.

$$p(x):=\sum _{i=1}^{n-1}(x-x_i{)}^2\in {\mathrm{\Pi }}_{2n-1}.$$

Hence this contradicts

$$\begin{array}{rcl}0& <& {\int }_a^b\omega (x)p(x)dx\\ & =& \sum _{i=1}^n{w}_{i}p(x_i)\\ & =& 0.\end{array}$$

For all $k=0,\dots ,n-1$,

$$\begin{array}{rcl}\sum _{i=1}^n{w}_i{p}_k(x_i)& =& {\int }_a^b\omega (x)p_k(x)dx\\ & =& (p_k,p_0)\\ & =& \{\begin{array}{ll}(p_0,p_0)& k=0\\ 0& 1\le k\le n-1\end{array}.\end{array}$$

Thus, $w_i$ satisfy $\text{(16)}$.

Next we will show $x_i$ are roots of $p_n$. Applying $\text{(17)}$ we will have

$$\begin{array}{rcl}0& =& (p_k,p_n)(\because \text{orthonormal})\\ & =& {\int }_a^b\omega (x)p_k(x)p_n(x)dx\\ \text{(20)}& & =& \sum _{i=1}^n{w}_i{p}_k(x_i)p_n(x)(\because \text{(17)}).\end{array}$$

In other word,

$$\begin{array}{r}A\left(\begin{array}{c}{w}_1{p}_n(x_1)\\ ⋮\\ w_n{p}_n(x_n)\end{array}\right)=0.\end{array}$$

By theorem 3.6.11, $A$ is nonsingular. Hence $c=0$, that is, for all $i=1,\dots ,n$,

$$w_i{p}_n(x_i)=0.$$

Since $w_i>0$, we have $p_n(x_i)=0$.

$◻$

Theorem 3.6.20

• ${\delta }_i,{\gamma }_i(i=1,\dots ,n)$,
  • defined in Theorem 3.6.3
• $p_0,\dots ,p_n$,
  • orthonomal polynomials
• $x_1,\dots ,x_n$,
  • the roots of $p_n$,

$$\begin{array}{r}{J}_j:=\left(\begin{array}{ccccc}{\delta }_1& {\gamma }_2& 0& & \\ {\gamma }_2& {\delta }_2& {\gamma }_3& & \\ 0& \cdot & \cdot & \cdot & \\ & & {\gamma }_j& {\delta }_{j-1}& {\gamma }_j\\ & & 0& {\gamma }_j& {\delta }_j\end{array}\right).\end{array}$$

Then $x_1,\dots ,x_n$ are the eigenvalues of the $J_n$.

proof

$$p_j(x):=\det (J_j-xI).$$

$◻$

Theorem 3.6.24

• $f\in C^{2n}[a,b]$,

$${\int }_a^b\omega (x)f(x)dx-\sum _{i=1}^n{w}_{i}f(x_i)=\frac{f^{(2n)}(\xi )}{(2n)!}(p_n,p_n).$$

proof

$◻$

Type of Gaussian quadrature

Gauss-Legendre

$$\begin{array}{}\text{(21)}& p_k(x):=\frac{k!}{(2k)1}\frac{d^k}{d{x}^k}(x-1{)}^k,k=0,1,\dots ,.\end{array}$$

• $\omega (x)=1$.
• $[a,b]:=[-1,1]$,

Gauss-Chebyshev

• $p_k$
  • Chebyshev polynomials
• $\omega (x):=(1-x^2{)}^{-1/2}$,
• $[a,b]:=[-1,1]$,

Gauss-Laguerre

• $p_k$
  • Laguerre polynomials
• $\omega (x):=e^{-x}$,
• $[a,b]:=[0,\mathrm{\infty }]$,

Gauss-Hermite

• $p_k$
  • Hermite polynomials
• $\omega (x):=e^{-x^2}$,
• $[a,b]:=[-\mathrm{\infty },\mathrm{\infty }]$,

$$\begin{array}{rcl}{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{1}{\sigma \sqrt{2\pi }}h(y)\exp (-\frac{(y-\mu {)}^2}{2{\sigma }^2})dy& =& {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{1}{\sigma \sqrt{2\pi }}h(\sqrt{2}\sigma x+\mu )\exp (-x^2)\sqrt{2}\sigma dy(x=(y-\mu )/\sqrt{2{\sigma }^2})\\ & =& {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\frac{1}{\sqrt{\pi }}h(\sqrt{2}\sigma x+\mu )\exp (-x^2)dy.\end{array}$$

Gauss-Jacobi

• $p_k$
  • Hermite polynomials
• $\omega (x):=(1-x{)}^{\alpha }(1+x{)}^{\beta }$,
• $[a,b]:=[-1,1]$,

Gauss-Radau

• In naive Gaussian quadrature, all abscissas are roots of the orthonomal polynomial. Gauss-Radau allows you to use preassigned nodes.
• A preassigned node is an endpoint of the interval, that is either $a$ or $b$.

Gauss-Lobatto

• In naive Gaussian quadrature, all abscissas are roots of the orthonomal polynomial. Gauss-Radau allows you to use preassigned nodes.
• preassigned nodes are endpoints of the interval, that is either $a$ and $b$.

Gauss-Kronrod

• In naive Gaussian quadrature, $x_i$ and $w_i$ need to be re-calcualted as $n$ increase. Gauss-Kronrod allows you to use the values calcualted in $n-1$.

Reference

• Stoer, Josef, and Roland Bulirsch. Introduction to numerical analysis. Vol. 12. Springer Science & Business Media, 2013.
• Hermite-Gauss Quadrature -- from Wolfram MathWorld
• Hermite polynomials - Wikipedia
• Gauss–Laguerre quadrature - Wikipedia