MATH566 Lesson 19: Numerical integration - polynomials

Overview

• Integration of polynomial approximant

• Integration operator: continuous and discrete cases, integration matrix

• Definite integration numerical formulas:

  • trapezoid

  • Simpson

  • Second Simpson

• Composite quadrature

Integration of polynomial approximants

• Consider $f:ℝ\to ℝ$, difficult to compute, and the problem of computation of

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

  The integration operator $I_{\mathrm{ab}}$ is linear $I_{\mathrm{ab}}(\alpha f+\beta g)=\alpha I_{\mathrm{ab}}\left(f\right)+\beta I_{\mathrm{ab}}\left(g\right)$

  ${\displaystyle \underset{a}{\overset{b}{\int }}[\alpha f\left(t\right)+\beta g\left(t\right)]dt=\alpha \underset{a}{\overset{b}{\int }}f\left(t\right)dt+\beta \underset{a}{\overset{b}{\int }}g\left(t\right)dt}$

• Approach to numerical integration is similar to numerical differentiation

  • Construct approximant of $f$, $g\cong f$

  • Approximate $I_{\mathrm{ab}}\left(f\right)$ by $I_{\mathrm{ab}}\left(g\right)$, $I_{\mathrm{ab}}\left(f\right)\cong I_{\mathrm{ab}}\left(g\right)$

• Easiest approximant to work with: $g$ polynomial, e.g., polynomial interpolant

  ${\displaystyle f\left(t\right)\cong g\left(t\right)=a_0+a_{1}t+\cdots +a_n{t}^n\Rightarrow I_{\mathrm{ab}}\left(f\right)=a_0(b-a)+\cdots +a_n\frac{b^{n+1}-a^{n+1}}{n+1}}$

Continuous and discrete integration operators

• $f:ℝ\to ℝ$, dificult to compute, only (partially) known through sample

  ${\displaystyle 𝒟=\{(x_i,f_i),i=0,\dots ,n\},f_i=f\left(x_i\right).}$

• Continuous-discrete operator analogy

  ${\displaystyle \begin{array}{ll}\mathrm{Continuous}& \mathrm{Discrete}\\ f:ℝ\to ℝ& 𝒟=\{(x_i,y_i=f\left(x_i\right)),i=0,\dots ,n\}\\ & 𝒙={\left[\begin{array}{lll}{x}_0& \dots & x_n\end{array}\right]}^T,𝒇={\left[\begin{array}{lll}{f}_0& \dots & f_n\end{array}\right]}^T\\ I_{\mathrm{ab}}:𝒞[a,b]\to ℝ& Q_{\mathrm{ab}}:{ℝ}^{n+1}\to ℝ\\ f\begin{array}{l}{I}_{\mathrm{ab}}\\ \to \\ \end{array}\underset{a}{\overset{b}{\int }}f\left(t\right)dt& 𝒇\begin{array}{l}{Q}_{\mathrm{ab}}\\ \to \\ \end{array}\sum _{i=0}^n{w}_i{f}_i\\ {\int }_a^{b}f\left(t\right)dt& {\sum }_{i=0}^n{w}_i{f}_i\cdot 1,e.g.,h=\frac{b-a}{n}\\ & \mathrm{Left}\mathrm{Darboux}:w_i=h,i=0,\dots ,n-1,w_n=0\\ & \mathrm{Right}\mathrm{Darboux}:w_0=0,w_i=h,i=0,\dots ,n\\ & \mathrm{Note}:\mathrm{multiple}\mathrm{possible}\mathrm{choices}\mathrm{for}{w}_i\end{array}}$

Integration (quadrature) matrix

• Computation of a definite integral yields a scalar from a function sample

  ${\displaystyle Q_{\mathrm{ab}}:{ℝ}^{n+1}\to ℝ,Q\left(𝒇\right)={𝒘}^T𝒇,𝒘\in {ℝ}^{n+1}}$

• Consider computation of the primitive

  ${\displaystyle F\left(x\right)=\underset{a}{\overset{x}{\int }}f\left(t\right)dt}$

  $f$ difficult to compute implies $F$ difficult to compute, known through sample

  ${\displaystyle {𝒟}^{\text{'}}=\{(x_i^{\text{'}},F_i^{\text{'}}=F\left(x_i\right)),i=0,1,\dots ,m\}}$

• Discrete analog, a linear mapping from ${ℝ}^{n+1}$ ($f$ sample) to ${ℝ}^{m+1}$ ($F$ sample)

  ${\displaystyle 𝑭={𝑲}_{𝒙,{𝒙}^{\text{'}}}𝒇}$ ${\displaystyle 𝑭=\left[\begin{array}{lll}{F}_0& \dots & F_m\end{array}\right],𝒇={\left[\begin{array}{lll}{f}_0& \dots & f_n\end{array}\right]}^T}$

Integration matrix components

• Lagrange form of polynomial interpolant of $𝒟=\{(x_i,f_i),i=0,\dots ,n\}$

  ${\displaystyle f\left(t\right)\cong p\left(t\right)=\sum _{j=0}^n{f}_j{\ell }_j\left(t\right),{\ell }_j\left(t\right)=\prod _{l=0}^n{}^{\text{'}}\frac{(t-x_l)}{(x_j-x_l)}}$

• Integral of polynomial interpolant

  ${\displaystyle F\left(x_i^{\text{'}}\right)\cong \underset{a}{\overset{x_i^{\text{'}}}{\int }}p\left(t\right)dt=\underset{a}{\overset{x_i^{\text{'}}}{\int }}\sum _{j=0}^n{f}_j{\ell }_j\left(t\right)dt=\sum _{j=0}^n(\underset{a}{\overset{x_i^{\text{'}}}{\int }}{\ell }_j\left(t\right)dt)f_j=\sum _{j=0}^n{k}_{ij}{f}_j}$

• In matrix form the above is

  ${\displaystyle 𝑭={𝑲}_{𝒙,{𝒙}^{\text{'}}}𝒇,{𝑲}_{𝒙,{𝒙}^{\text{'}}}\in {ℝ}^{(m+1)\times (n+1)}}$

• Components of ${𝑲}_{𝒙,{𝒙}^{\text{'}}}$

  ${\displaystyle k_{ij}=(\underset{a}{\overset{x_i^{\text{'}}}{\int }}{\ell }_j\left(t\right)dt)}$

Definite integration quadrature formulas

• To approximate $I_{\mathrm{ab}}\left(f\right)={\int }_a^{b}f\left(t\right)dt$, evaluate $Q_{\mathrm{ab}}\left(𝒇\right)\cong I_{\mathrm{ab}}\left(f\right)$

  ${\displaystyle Q_{\mathrm{ab}}:{ℝ}^{n+1}\to ℝ,Q_{\mathrm{ab}}\left(𝒇\right)={𝒘}^T𝒇,𝒘\in {ℝ}^{n+1}}$

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

• Construct polynomial interpolant $f\left(t\right)\cong p\left(t\right)={\sum }_{i=0}^n{f}_i{\ell }_i\left(t\right)$

  ${\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}$

• Trapezoid formula: $n=1$, $p$ linear, $b=x_1=x_0+h=a+h$,

  ${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong (\underset{x_0}{\overset{x_1}{\int }}\frac{t-x_1}{x_0-x_1}dt)f_0+(\underset{x_0}{\overset{x_1}{\int }}\frac{t-x_0}{x_1-x_0}dt)f_1\Rightarrow }$ ${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong Q_{\mathrm{ab}}\left(f\right)=\frac{h}{2}{f}_0+\frac{h}{2}{f}_1=\frac{b-a}{2}(f_0+f_1)={𝒘}^T𝒇,𝒘=\left[\begin{array}{ll}h/2& h/2\end{array}\right].}$

Definite integration quadrature formulas

• Simpson formula: $n=2$, $p$ quadratic, $x_0=a,x_1=a+h,x_2=a+2h=b$,

  ${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong Q_{\mathrm{ab}}\left(𝒇\right)=\sum _{i=0}^n(\underset{a}{\overset{b}{\int }}{\ell }_i\left(t\right)dt)f_i}$ ${\displaystyle w_0=\underset{a}{\overset{b}{\int }}{\ell }_0\left(t\right)dt=\underset{x_0}{\overset{x_2}{\int }}\frac{(t-x_1)(t-x_2)}{(x_0-x_1)(x_0-x_2)}dt=\frac{h}{3}}$ ${\displaystyle w_1=\underset{a}{\overset{b}{\int }}{\ell }_1\left(t\right)dt=\underset{x_0}{\overset{x_2}{\int }}\frac{(t-x_0)(t-x_2)}{(x_1-x_0)(x_1-x_2)}dt=\frac{4h}{3}}$ ${\displaystyle w_2=\underset{a}{\overset{b}{\int }}{\ell }_2\left(t\right)dt=\underset{x_0}{\overset{x_2}{\int }}\frac{(t-x_0)(t-x_1)}{(x_2-x_0)(x_2-x_1)}dt=\frac{h}{3}}$ ${\displaystyle }$ ${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong Q_{\mathrm{ab}}\left(𝒇\right)=\frac{h}{3}(f_0+4f_1+f_2)=\frac{b-a}{6}(f_0+4f_1+f_2)}$ ${\displaystyle {𝒘}^T=\frac{h}{3}\left[\begin{array}{lll}1& 4& 1\end{array}\right]}$

Definite integration quadrature formulas

• Second Simpson formula: $n=3$, $p$ cubic, $x_i=a+ih,h=(b-a)/3$

  ${\displaystyle \underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong Q_{\mathrm{ab}}\left(𝒇\right)=\sum _{i=0}^n(\underset{a}{\overset{b}{\int }}{\ell }_i\left(t\right)dt)f_i}$ ${\displaystyle w_0=\underset{x_0}{\overset{x_3}{\int }}{\ell }_0\left(t\right)dt=\underset{x_0}{\overset{x_3}{\int }}\frac{(t-x_1)(t-x_2)(t-x_3)}{(x_0-x_1)(x_0-x_2)(x_0-x_3)}dt=\frac{3h}{8}}$
  ${\displaystyle w_1=\underset{a}{\overset{b}{\int }}{\ell }_1\left(t\right)dt=\underset{x_0}{\overset{x_3}{\int }}\frac{(t-x_0)(t-x_2)(t-x_3)}{(x_1-x_0)(x_1-x_2)(x_1-x_3)}dt=\frac{9h}{8}}$ ${\displaystyle w_2=\underset{a}{\overset{b}{\int }}{\ell }_2\left(t\right)dt=\underset{x_0}{\overset{x_3}{\int }}\frac{(t-x_0)(t-x_1)(t-x_3)}{(x_2-x_0)(x_2-x_1)(x_2-x_3)}dt=\frac{9h}{8}}$ ${\displaystyle w_3=\underset{a}{\overset{b}{\int }}{\ell }_3\left(t\right)dt=\underset{x_0}{\overset{x_3}{\int }}\frac{(t-x_0)(t-x_1)(t-x_2)}{(x_3-x_0)(x_3-x_1)(x_3-x_2)}dt=\frac{3h}{8}}$ ${\displaystyle }$ ${\displaystyle Q_{\mathrm{ab}}\left(𝒇\right)=\frac{h}{8}(3f_0+9f_1+9f_2+3f_3)=\frac{b-a}{8}(f_0+3f_1+3f_2+f_2)}$

Error analysis

• Previous quadrature formulas based upon $k$-degree polynomial interpolation

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)=\underset{a}{\overset{b}{\int }}f\left(t\right)dt\cong Q_{\mathrm{ab}}^{\left(k\right)}\left(f\right)=(b-a)\sum _{i=0}^k{w}_i^{\left(k\right)}{f}_i,f_i=f(a+ih),h=\frac{b-a}{k+1}}$ ${\displaystyle \begin{array}{llll}k& 1& 2& 3\\ \mathrm{Name}& \mathrm{Trapezoid}& \mathrm{Simpson}& 2^{\mathrm{nd}}\mathrm{Simpson}\\ 𝒘& \frac{1}{2}\left[\begin{array}{l}1\\ 1\end{array}\right]& \frac{1}{6}\left[\begin{array}{l}1\\ 4\\ 1\end{array}\right]& \frac{1}{8}\left[\begin{array}{l}1\\ 3\\ 3\\ 1\end{array}\right]\end{array}}$

• Numerical integration error is evaluated by integral of interpolation error

  ${\displaystyle e=I_{\mathrm{ab}}\left(f\right)-Q_{\mathrm{ab}}^{\left(k\right)}\left(f\right)=\underset{a}{\overset{b}{\int }}[f\left(t\right)-p_k\left(t\right)]dt=\underset{a}{\overset{b}{\int }}\frac{f^{(k+1)}\left({\xi }_t\right)}{(k+1)!}w\left(t\right)dt}$ ${\displaystyle w\left(t\right)=\prod _{i=0}^k(t-x_i),x_i=a+ih}$

Common single-interval error estimates

• Trapezoid

  ${\displaystyle e=|\frac{1}{2}\underset{a}{\overset{b}{\int }}{f}^{{}^{\text{'}\text{'}}}\left({\xi }_t\right)(t-a)(t-b)dt|⩽{||f^{\text{'}\text{'}}||}_{\infty }\frac{1}{2}|\underset{-h/2}{\overset{h/2}{\int }}(y^2-{\left(\frac{h}{2}\right)}^2)dy|\Rightarrow }$ ${\displaystyle e⩽{||f^{\text{'}\text{'}}||}_{\infty }\underset{0}{\overset{h/2}{\int }}({\left(\frac{h}{2}\right)}^2-y^2)dy={||f^{\text{'}\text{'}}||}_{\infty }\left(\frac{h^3}{12}\right)}$

• Simpson

  • Try with $k=2$ to obtain $e⩽{||f^{\left(3\right)}||}_{\infty }\frac{1}{6}{\int }_a^b(t-a)(t-\frac{a+b}{2})(t-b)dt$

    ${\displaystyle \underset{a}{\overset{b}{\int }}(t-a)(t-\frac{a+b}{2})(t-b)dt=\underset{-h}{\overset{h}{\int }}(y-h)y(y+h)dt=0}$

  • Terms of order $𝒪\left(h^4\right)$ cancel out!

  • Peano kernel theorem is used to obtain

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

Composite quadrature

• Increased accuracy arises from more sample points $𝒟=\{(x_i,f_i),i=0,\dots ,m\}$

  • Do not use high degree polynomial interpolant (Runge phenomenon)

  • Apply a lower-order formula repeatedly over subintervals $\Rightarrow$ “composite”

• Composite trapezoid formula

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)=\underset{a}{\overset{b}{\int }}f\left(t\right)dt=\sum _{i=1}^m\underset{x_{i-1}}{\overset{x_i}{\int }}f\left(t\right)dt}$

  Assume equidistant sampling $x_{i+1}-x_i=h=(b-a)/m$

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)\cong Q_{\mathrm{ab}}\left(f\right)=\frac{h}{2}\sum _{i=1}^m(f_{i-1}+f_i)=\frac{h}{2}(f_0+f_m)+h\sum _{i=1}^{m-1}{f}_i}$

• Composite trapezoid error, over $[x_{i-1},x_i]$, error is $e_i⩽{||f^{\text{'}\text{'}}||}_{\infty }{h}^3/12$

  ${\displaystyle e⩽\sum _{i=1}^m{e}_i⩽{||f^{\text{'}\text{'}}||}_{\infty }m{h}^3/12={||f^{\text{'}\text{'}}||}_{\infty }\frac{b-a}{12}{h}^2}$

  Composite trapezoid is second-order accurate.

Compiste Simpson

• Assume $𝒟=\{(x_i,f_i),i=0,1,\dots ,m\}$, $m=2p$

• Composite Simpson formula

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)=\underset{a}{\overset{b}{\int }}f\left(t\right)dt=\sum _{i=1}^p\underset{x_{2i-2}}{\overset{x_{2i}}{\int }}f\left(t\right)dt}$

  Assume equidistant sampling $h=(b-a)/m=(b-a)/\left(2p\right),x_{2i}-x_{2i-2}=2h$

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)\cong Q_{\mathrm{ab}}\left(f\right)=\frac{h}{3}\sum _{i=1}^p(f_{2i-2}+4f_{2i-1}+f_{2i})}$ ${\displaystyle }$

• Composite Simpson error, over $[x_{2i-2},x_{2i}]$, error is $e_i⩽{||f^{\left(4\right)}||}_{\infty }{h}^5/90$

  ${\displaystyle e⩽\sum _{i=1}^m{e}_i⩽{||f^{\left(4\right)}||}_{\infty }{h}^5/90={||f^{\left(4\right)}||}_{\infty }\frac{b-a}{180}{h}^4}$

  Composite Simpson is fourth-order accurate.