MATH566 Lesson 21: Romberg quadrature

Overview

• Sequence convergence acceleration

• Aitken acceleration

• Composite trapezoid + Aitken acceleration = Romberg quadrature

Sequence convergence acceleration

• Recall definitions that characterize sequence convergence

  Definition. ${\left\{x_n\right\}}_{n\in \mathbb{N}}$ converges to $x$ with rate $r\in (0,1)$ and order $p$ if

  ${\displaystyle {lim}_{n\to \infty }\frac{|x_{n+1}-x|}{{|x_n-x|}^p}=r.}$

  (1)

• If $p$ is known, it can be used to construct a faster converging sequence $\left\{y_n\right\}$

• Write $x_n-x\cong r{(x_{n-1}-x)}^p,x_{n-1}-x\cong r{(x_{n-2}-x)}^p$, and take ratio

  ${\displaystyle \frac{x_n-\bar{x}}{x_{n-1}-\bar{x}}={\left(\frac{x_{n-1}-\bar{x}}{x_{n-2}-\bar{x}}\right)}^p.}$

  (2)

• If $p,x_{n-2},x_{n-1},x_n$ are known, then the above can be solved for $\bar{x}$

• The solution $\bar{x}$ is taken as the $n^{\mathrm{th}}$ term $y_n=\bar{x}$ in a faster converging sequence

• Since it is based upon assumed behavior, this is called an extrapolation, understood as “going outside realm of known data”

• Other extrapolation example: Let $p_m\left(t\right)$ interpolate data $𝒟=\{(ih,y_i),i=0,1,\dots ,m\}$. $p_m(jh)$ for $j\in [0,m]$ is interpolation, for $j>m$ is an extrapolation.

Aitken acceleration

• Extrapolation for $p=1$ is known as Aitken acceleration

  ${\displaystyle \frac{x_n-\bar{x}}{x_{n-1}-\bar{x}}=\frac{x_{n-1}-\bar{x}}{x_{n-2}-\bar{x}}}$

• Obtain

  ${\displaystyle x_n{x}_{n-2}-(x_n+x_{n-2})\bar{x}=x_{n-1}^2-2x_{n-1}\bar{x}\Rightarrow \bar{x}=\frac{x_n{x}_{n-2}-x_{n-1}^2}{x_n-2x_{n-1}+x_{n-2}}.}$

• Set $y_n=\bar{x}$, rewritten as a correction to the $x_n$ term

  ${\displaystyle y_n=x_n-\frac{{(x_n-x_{n-1})}^2}{x_n-2x_{n-1}+x_{n-2}}}$

• Example: L21.jl shows acceleration of forward finite difference approximation

  ${\displaystyle f_0^{\text{'}}\cong \frac{f(x_0+h_n)-f\left(x_0\right)}{h_n},h_n=2^{-n}}$

Composite trapezoid

• Consider problem of computing $I_{\mathrm{ab}}\left(f\right)={\int }_a^{b}f\left(t\right)dt$ by composite trapezoid with uniform partition of interval $[a,b]$, $x_i=x_0+i{h}_n$, $h_n=(b-a)/n$

  ${\displaystyle I_{\mathrm{ab}}\left(f\right)\cong Q_{\mathrm{ab}}\left(f\right)=T\left(n\right)=(h_n/2)[f_0+2f_1+2f_2+\cdots +2f_{n-1}+f_n]}$

• Introduce notation ${\sum }^{\text{'}\text{'}}$to indicate that first, last terms in sum are halved

  ${\displaystyle \sum _{i=0}^n{}^{\text{'}\text{'}}{f}_i=\frac{1}{2}{f}_0+f_1+f_2+\cdots +f_{n-1}+\frac{1}{2}{f}_n,T\left(n\right)=h_n\sum _{i=0}^n{}^{\text{'}\text{'}}{f}_i}$

• Recall: computational complexity of numerical quadrature given by number $f$ evaluations. Consider $T\left(n\right)$,$T\left(2n\right)$, and let $f_i^{\left(n\right)}=f(i{h}_n)$

  ${\displaystyle \begin{array}{llllllllll}T\left(n\right)& =& h_n\cdot & (.& f_0^{\left(n\right)}+& & f_1^{\left(n\right)}+\cdots & +& f_n^{\left(n\right)}& .)\\ T\left(2n\right)& =& h_{2n}\cdot & (.& f_0^{(n+1)}+& f_1^{(n+1)}+& f_2^{(n+1)}+\dots & +& f_{2n}^{\left(2n\right)}& .)\end{array}}$

• Many common points: exploit this to reduce $f$ evaluations

Recursive composite trapezoid formula

• Consider $T\left(n\right)$ to have been computed, e.g.,

  ${\displaystyle T\left(1\right)=\frac{b-a}{2}(f\left(a\right)+f\left(b\right))}$

• Let $h=(b-a)/\left(2n\right)$, express $T\left(2n\right)$ in terms of previously computed $T\left(n\right)$

  ${\displaystyle T\left(2n\right)=\frac{1}{2}T\left(n\right)+h\sum _{i=1}^n{f}_{2i-1},f_j=f(jh),x_j=a+jh,j=0,\dots ,2n.}$

• The above allows construction of sequence ${\left\{x_n\right\}}_{n\in \mathbb{N}}$, $x_n=T\left(n\right)$.

• Recall composite trapezoid error formula

  ${\displaystyle e_n=|I_{\mathrm{ab}}\left(f\right)-T\left(n\right)|⩽{||f^{\text{'}\text{'}}||}_{\infty }\frac{b-a}{12}{h}_n^2}$

  of second-order accuracy

Romberg integration

• Idea #1: apply extrapolation acceleration to $\left\{x_n\right\}$ sequence

  ${\displaystyle T\left(n\right)-I_{\mathrm{ab}}\cong {||f^{\text{'}\text{'}}||}_{\infty }\frac{b-a}{12}{h}_n^2,T\left(2n\right)-I_{\mathrm{ab}}\cong {||f^{\text{'}\text{'}}||}_{\infty }\frac{b-a}{12}{h}_{2n}^2\Rightarrow }$ ${\displaystyle \frac{T_n-I_{\mathrm{ab}}}{T_{2n}-I_{\mathrm{ab}}}=4\Rightarrow I_{\mathrm{ab}}=T\left(2n\right)+\frac{1}{4-1}[T\left(2n\right)-T\left(n\right)]}$

• Take $y_n=I_{\mathrm{ab}}=T\left(2n\right)+\frac{1}{4-1}[T\left(2n\right)-T\left(n\right)]$, which should now be faster-convergent, (third-order) sequence

• Idea #2: repeatedly apply extrapolation acceleration to $\left\{y_n\right\}$ sequence

• Introduce $R(n,q)$, $n$ = number of subintervals, $q$ = number of times extrapolation acceleration has been applied

• Romberg relation

  ${\displaystyle R(n,q)=R(n,q-1)+\frac{1}{4^q-1}[R(n,q-1)-R(n-1,q-1)]}$

Romberg table

• Organize computations as a table

  ${\displaystyle \begin{array}{llllll}\mathrm{Resolution}& \mathrm{Extrapolation}\mathrm{steps}\to & & & & \\ \downarrow & R(0,0)& & & & \\ & R(1,0)& R(1,1)& & & \\ & R(2,0)& R(2,1)& R(2,2)& & \\ & ⋮& ⋮& ⋮& \ddots & \\ & R(n,0)& R(n,1)& R(n,2)& \dots & R(n,q)\end{array}}$

• Choose $n$ high enough to ensure resolution, i.e., function variation is captured

• Choose moderate $q$, $q⩽4$, to avoid loss of precision due to floating point subtraction