MATH566 Lesson 22: Gauss quadrature

Overview

• Weighted quadrature - method of moments

• Optimal sampling - Gauss quadrature

• Common Gauss quadrature methods:

  • Gauss-Legendre ${\int }_{-1}^{1}f\left(t\right)dt$

  • Gauss-Chebyshev ${\int }_{-1}^1\frac{f\left(t\right)}{\sqrt{1-t^2}}dt$

  • Gauss-Laguerre ${\int }_0^{\infty }{e}^{-t}f\left(t\right)dt$

  • Gauss-Legendre ${\int }_{-\infty }^{\infty }{e}^{-t^2}f\left(t\right)dt$

Method of moments: review of simplest case

• $f:ℝ\to ℝ$ known through data set ($f$ sample) $𝒟=\{(x_i,f_i),i=0,1,\dots ,n\}$

  ${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong \sum _{i=0}^n(\underset{a}{\overset{b}{\int }}{\ell }_i\left(t\right)dt)f_i=\sum _{i=0}^n{w}_i{f}_i}$

• Set a simple, predefined integration domain, e.g., $[0,1]$, $x_i=ih$, $h=1/n$

• Monomial basis set: $ℳ=\{1,t,t^2,\dots \}$ conditions

  ${\displaystyle \begin{array}{llllll}f\left(t\right)=1:& \underset{0}{\overset{1}{\int }}1dt& =& 1& =& \sum _{i=0}^n{w}_i\\ f\left(t\right)=t:& \underset{0}{\overset{1}{\int }}tdt& =& \frac{1}{2}& =& h\sum _{i=0}^n{w}_{i}i\\ f\left(t\right)=t^2:& \underset{0}{\overset{1}{\int }}{t}^{2}dt& =& \frac{1}{3}& =& h\sum _{i=0}^n{w}_i{i}^2\\ \dots .& & & & & \end{array}}$

• Solve above system to find weights $w_i$, obtain formula with error $e\propto {||f^{\left(n\right)}||}_{\infty }$

• What if $f^{\left(n\right)}$ is not defined (infinite) at some points in the integration domain?

Integrals with singular integrands

• Consider integrals of form

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)=\underset{a}{\overset{b}{\int }}\omega \left(t\right)f\left(t\right)dt}$

• Typically:

  • $f$ is smooth, e.g., $f\in C^{\infty }[a,b]$ (all derivatives exist, are finite)

  • $w$ captures some singular behavior of the integral

• Examples:

  • Chebyshev weight

    ${\displaystyle T\left(f\right)=\underset{-1}{\overset{1}{\int }}\frac{f\left(t\right)}{\sqrt{1-t^2}}dt}$

  • Laguerre weight

    ${\displaystyle L\left(f\right)=\underset{0}{\overset{\infty }{\int }}{e}^{-t}f\left(t\right)dt}$

  • Hermite weight

    ${\displaystyle H\left(f\right)=\underset{-\infty }{\overset{\infty }{\int }}{e}^{-t^2}f\left(t\right)dt}$

Method of moments for weighted integrals

• As in the $\omega \left(t\right)=1$ case, choose $f$ evaluation points, $x_i=a+ih$,$h=(b-a)n$

• Choose a basis set, e.g., monomials $ℳ=\{1,t,t^2,\dots \}$ and impose exact result when using exact, analytical weight function $\omega \left(t\right)$

• Example for Chebyshev, $\omega \left(t\right)={(1-t^2)}^{-1/2}$, $x_i=ih-1$, $h=2/n$, $i=0,\dots ,n$

  ${\displaystyle \begin{array}{llllllll}f\left(t\right)=1:& \underset{-1}{\overset{1}{\int }}\omega \left(t\right)1dt& =& \underset{-1}{\overset{1}{\int }}\frac{1}{\sqrt{1-t^2}}dt& =& \pi & =& \sum _{i=0}^n{w}_i\\ f\left(t\right)=t:& \underset{-1}{\overset{1}{\int }}\omega \left(t\right)tdt& =& \underset{-1}{\overset{1}{\int }}\frac{t}{\sqrt{1-t^2}}dt& =& 0& =& \sum _{i=0}^n{w}_i(ih-1)\\ f\left(t\right)=t^2:& \underset{-1}{\overset{1}{\int }}\omega \left(t\right)t^{2}dt& =& \underset{-1}{\overset{1}{\int }}\frac{t^2}{\sqrt{1-t^2}}dt& =& \frac{\pi }{2}& =& \sum _{i=0}^n{w}_i{(ih-1)}^2\\ \dots .& & & & & & & \end{array}}$

Optimal sampling

• Up to now the sampling points $x_i$ were arbitrarily chosen

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)=\underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong Q_{\mathrm{ab}}\left(f\right)=\sum _{i=0}^n{w}_{i}f\left(x_i\right)}$

• One could obtain more accurate method by optimal choice of $x_i$, more of the moment equations could be satisfied, $2(n+1)$ instead of only $n+1$

• Difficult to do within method of moments since a non-linear system results

${\displaystyle \begin{array}{llllll}f\left(t\right)=1:& \underset{0}{\overset{1}{\int }}1dt& =& 1& =& \sum _{i=0}^n{w}_i\\ f\left(t\right)=t:& \underset{0}{\overset{1}{\int }}tdt& =& \frac{1}{2}& =& \sum _{i=0}^n{w}_i{x}_i\\ f\left(t\right)=t^2:& \underset{0}{\overset{1}{\int }}{t}^{2}dt& =& \frac{1}{3}& =& \sum _{i=0}^n{w}_i{x}_i^2\\ \dots .& & & & & \end{array}}$

Gauss quadrature

• Solving the nonlinear system can be avoided by use of orthogonal polynomials

• Consider

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)=\underset{a}{\overset{b}{\int }}\omega \left(t\right)f\left(t\right)dt\cong Q_{\mathrm{ab}}\left(f\right)=\sum _{i=0}^n{w}_i{f}_i,f_i=f\left(x_i\right)}$

• In method of moments there are $2(n+1)$ parameters, $n+1$ weights, $n+1$ $x_i$'s

• Can impose exact quadrature up to degree $2n+1$

• Consider $p_{2n+1}\left(t\right)$ a polynomial of degree $2n+1$

  ${\displaystyle I_{\mathrm{ab}}\left(p\right)=\underset{a}{\overset{b}{\int }}\omega \left(t\right)p\left(t\right)dt}$

• Introduce scalar product

  ${\displaystyle {(f,g)}_{\omega }=\underset{a}{\overset{b}{\int }}\omega \left(t\right)f\left(t\right)g\left(t\right)dt}$

  with $\omega \left(t\right)$ assumed to satisfy scalar product properties ($t\in [a,b]\Rightarrow \omega \left(t\right)⩾0$)

Gauss quadrature: polynomial division

• Introduce ${\phi }_k\left(t\right)$ polynomials orthogonal w.r.t. ${(f,g)}_{\omega }$ scalar product

  ${\displaystyle {({\phi }_j,{\phi }_k)}_{\omega }={\delta }_{jk}}$

• ${\phi }_k$ polynomials can be found by Gram-Schmidt algorithm applied to $\{1,t,\dots \}$

• Divide $p_{2n+1}\left(t\right)$ by ${\phi }_{n+1}\left(t\right)$

  ${\displaystyle p_{2n+1}\left(t\right)=q_n\left(t\right){\phi }_{n+1}\left(t\right)+r_n\left(t\right)}$

• Example $p\left(t\right)=3t^3-2t^2+t+1$ divided by ${\phi }_1\left(t\right)=t^2+t+1$, $p_3\left(t\right)=q_1\left(t\right){\phi }_2\left(t\right)+r_1\left(t\right)$, $q_1\left(t\right)=3t-5$, $r_1\left(t\right)=3t+6$

  ${\displaystyle \begin{array}{lllll}& 3t& -5& & \\ t^2+t+1& 3t^3& -2t^2& +t& +1\\ & 3t^3& +3t^2& +3t& \\ & & -5t^2& -2t& +1\\ & & -5t^2& -5t& -5\\ & & & 3t& +6\end{array}}$

Gauss quadrature: basic idea

• Integrate $p_{2n+1}$

  ${\displaystyle \underset{a}{\overset{b}{\int }}\omega \left(t\right)p_{2n+1}\left(t\right)dt=\underset{a}{\overset{b}{\int }}\omega \left(t\right)q_n\left(t\right){\phi }_{n+1}\left(t\right)dt+\underset{a}{\overset{b}{\int }}\omega \left(t\right)r_n\left(t\right)dt}$

• However, ${\phi }_{n+1}$ orthogonal w.r.t. to all polynomials of degree at most $n\Rightarrow$

  ${\displaystyle \underset{a}{\overset{b}{\int }}\omega \left(t\right)q_n\left(t\right){\phi }_{n+1}\left(t\right)dt=0}$

• Obtain

  ${\displaystyle \underset{a}{\overset{b}{\int }}\omega \left(t\right)p_{2n+1}\left(t\right)dt=\underset{a}{\overset{b}{\int }}\omega \left(t\right)r_n\left(t\right)dt=\sum _{i=0}^n{w}_i{r}_n\left(x_i\right)}$

  implying that a formula that is exact for polynomials up to degree $n$ is also exact for $p_{2n+1}$ of degree $2n+1$.

• Gauss observation: the only values of $p_{2n+1}$ that arise are $p_{2n+1}\left(x_i\right)$. Can $x_i$ be chosen such that $r_n\left(x_i\right)=p_{2n+1}\left(x_i\right)=q_n\left(x_i\right){\phi }_{n+1}\left(x_i\right)+r_n\left(x_i\right)$? Yes!

• Choose evaluation point $x_i$ such that ${\phi }_{n+1}\left(x_i\right)=0$, roots of ${\phi }_{n+1}$

Gauss-Legendre quadrature

• $I\left(f\right)={\int }_{-1}^{1}f\left(t\right)dt$, weight $\omega \left(t\right)=1$, scalar product $(f,g)={\int }_{-1}^{1}f\left(t\right)g\left(t\right)dt$

• Orthogonal family of polynomials are the Legendre polynomials

  ${\displaystyle \{1,t,\frac{1}{2}(3t^2-1),\frac{1}{2}(5t^3-3t),\frac{1}{8}(35t^4-30t^2+3),\dots \}}$

• Gauss-Legendre quadrature formulas

  • GL2 $n+1=2\Rightarrow {\phi }_2\left(t\right)=\frac{1}{2}(3t^2-1)$ with roots $x_0=-\frac{1}{\sqrt{3}}$, $x_1=\frac{1}{\sqrt{3}}.$

    ${\displaystyle I\left(f\right)=\underset{-1}{\overset{1}{\int }}f\left(t\right)dt\cong Q\left(f\right)=f\left(x_0\right)+f\left(x_1\right),w_0=w_1=1,}$

    exact up to cubics $I\left(p_3\right)=Q\left(p_3\right)$

  • GL3 $n+1=3\Rightarrow {\phi }_3\left(t\right)=\frac{1}{2}(5t^3-3t)$ with roots $x_{0,2}=-\frac{\sqrt{3}}{\sqrt{35}},x_1=0$,

    ${\displaystyle \underset{-1}{\overset{1}{\int }}f\left(t\right)dt\cong Q\left(f\right)=\frac{5}{9}f\left(x_0\right)+\frac{8}{9}f\left(x_1\right)+\frac{5}{9}f\left(x_2\right),w_0=w_2=\frac{5}{9},w_1=\frac{8}{9}}$

    exact up to quintics $I\left(p_5\right)=Q\left(p_5\right)$

Gauss-Chebyshev quadrature

• Gauss-Chebyshev sample points given by roots of $T_{n+1}\left(\theta \right)=\cos (n\theta )$, $x=\cos \theta$

  ${\displaystyle x_i=\cos \frac{(2i+1)\pi }{2(n+1)},i=0,1,2,\dots ,n}$

• Gauss-Chebyshev weights are especially simple, they're all equal!

  ${\displaystyle w_i=\frac{\pi }{n+1}}$ ${\displaystyle \underset{-1}{\overset{1}{\int }}\frac{f\left(t\right)}{\sqrt{1-t^2}}dt\cong Q\left(f\right)=\sum _{i=0}^n{w}_i{f}_i,f_i=f\left(x_i\right)}$

Gauss-Laguerre quadrature

• Relevant for Laplace transforms

  ${\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-t}f\left(t\right)dt\cong Q\left(f\right)=\sum _{i=0}^n{w}_i{f}_i,f_i=f\left(x_i\right)}$

• Gauss-Laguerre 2, exact up to cubics

  ${\displaystyle x_0=0.585786,w_0=0.853853,x_1=3.41421,w_1=0.146447}$

• Gauss-Laguerre 3, exact up to quintics

  ${\displaystyle x_0=0.415775,w_0=0.711093,x_1=2.29428,w_1=0.278518}$ ${\displaystyle x_2=6.28995,w_2=0.0103893}$

Gauss-Hermite quadrature

• Relevant for Gauss distributions, diffusion equations

  ${\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{e}^{-t^2}f\left(t\right)dt\cong Q\left(f\right)=\sum _{i=0}^n{w}_i{f}_i,f_i=f\left(x_i\right)}$

• Gauss-Hermite 2, exact up to cubics

  ${\displaystyle x_{0,1}=\pm 0.707107,w_{0,1}=0.886227}$

• Gauss-Hermite 3, exact up to quintics

  ${\displaystyle x_{0,2}=\pm 1.22474,w_{0,1}=0.0813128,x_1=0,w_1=1.18164}$

• Wiener transform $F\left(x\right)$ of $f\left(y\right)$ is

  ${\displaystyle F\left(x\right)=\frac{1}{\sqrt{4\pi }}\underset{-\infty }{\overset{\infty }{\int }}{e}^{-{(x-y)}^2/4}f\left(y\right)dy}$

  specifies the temperature after one unit of time of a bar whose initial temperature was $f\left(y\right)$.