Integration by Richardson Extrapolation

The Richardson extrapolation can also be applied to improve numerical integration (with resulting method called Romberg Method). Recall the composite trapezoidal rule in Eq. (71):

$$I[f] = \int_a^b f(x)\,dx \approx h\frac{f(a)+f(b)}{2}+h\,\sum_{k=1}^{n-1}f(a+kh)$$

$$= h\frac{f(a)-f(b)}{2}+h\,\sum_{k=1}^nf(a+kh)$$ (137)

where $h=(b-a)/n$. In theory, taking the limit $n\rightarrow\infty$ we get $h=(b-a)/n\rightarrow 0$, $I_0(h)$ becomes the true integral $I[f]$:

$$\lim\limits_{h\rightarrow 0}I_0(h)=\int_a^bf(x)\,dx=I[f]$$ (138)

Of course this approach is not practical computationally.

However, we can apply the Richardson extrapolation to increase the order of the error terms and improve the accuracy. Specifically, we introduce another variable $t$ such that $x=a+t\,h$, $dx=h\,dt,\;\;x_k=a+k\,h,\;\;b=a+nh$. We further define $\phi(t)=f(a+t\,h)$ and get

$$\phi^{(k)}(t)=\frac{d^k}{dt^k}\phi(t) =\frac{d^k}{dt^k}f(a+th)=h^kf^{(k)}(x)$$ (139)

Now $f(a)$, $f(a+k\,h)$, and $f(b)=f(a+nh)$ become $\phi(0)$, $\phi(k)$, and $\phi(n)$, respectively, and the composite trapezoidal method above becomes

$$I[f] = \int_a^b f(x)\,dx=h \int_0^n \phi(t)\,dt\approx h\frac{\phi(0)-\phi(n)}{2}+h\sum_{k=1}^n\phi(k)$$ (140)

Comparing this to the Euler-Maclaurin Formula:

$$\int_0^n \phi(t)\,dt=\sum_{k=1}^n\phi(k)+\frac{\phi(0)-\phi(n)}{2} -\sum_{k=1}^\infty\frac{B_{2k}}{(2k)!}[\phi^{(2k-1)}(n)-\phi^{(2k-1)}(0)]$$ (141)

we get

$$I[f] = \int_a^b f(x)\,dx=h \int_0^n \phi(t)\,dt$$

$$= h\frac{\phi(0)-\phi(n)}{2}+h\sum_{k=1}^n\phi(k) -h\sum_{k=1}^\infty\frac{B_{2k}}{(2k)!}[\phi^{(2k-1)}(n)-\phi^{(2k-1)}(0)]$$

$$= h\left[\frac{f(a)-f(b)}{2}+\sum_{k=1}^nf(a+kh) \right] -h^{2k}\sum_{k=1}^\infty\frac{B_{2k}}{(2k)!}[f^{(2k-1)}(b)-f^{(2k-1)}(a)]$$

$$= I_0(h)+\sum_{k=1}^\infty C_k h^{2k}=I_0(h)+C_1 h^2+C_4 h^4+\cdots =I_0(h)+O(h^2)$$ (142)

where $I_0(h)$ denotes the method of the composite trapezoidal rule

$$I_0(h)=h\frac{f(a)-f(b)}{2}+h\sum_{k=1}^nf(a+kh)$$ (143)

with a step size $h=(b-a)/n$ and an error term $O(h^2)$, and

$$C_k=-\frac{B_{2k}}{(2k)!}[f^{(2k-1)}(b)-f^{(2k-1)}(a)], \;\;\;\;\;(k=1,\,2,\,3,\cdots)$$ (144)

Now we can use the Richardson extrapolation to improve the accuracy of this method $I_0(h)$, by replacing $h$ in Eq. (142) by $h/2$:

$$I[f]=I(h/2)+C_1 \left(\frac{h}{2}\right)^2 +C_2\left(\frac{h}{2}\right)^4+O(h^6)$$ (145)

and then Subtract Eq. (98) from 4 times this equation to get

$$I[f]=\frac{4 I(h/2)-I_n(h)}{3}+O(h^4)=I_1(h/2)+O(h^4)$$ (146)

where

$$I_1(h/2)=\frac{4 I(h/2)-I(h)}{3}$$ (147)

is a new method with a step size $h/2$ and an error term $O(h^4)$.

The accuracy can be further improved by the same process. Again replacing $h$ in $I_1(h/2)$ by $h/2$ and repeating the same steps we get

$$I[f]=\frac{16 I_1(h/4)-I_1(h/2)}{15}+O(h^6)=I_2(h/4)+O(h^6)$$ (148)

where we defined a new method:

$$I_2(h/4)=\frac{16 I_1(h/4)-I_1(h/2)}{15}$$ (149)

with a step size $h/4$ and an error term $O(h^6)$.

Example:

Find the integral of the Gaussian function with zero mean and unit standard deviation:

$$I[f]=\int_a^b f(x)\,dx=\frac{1}{\sqrt{2\pi}}\int_0^3 e^{-x^2/2}\;dx$$ (150)

Initially, we set $m=0$, $h_0=(b-a)/2^m=b-a=3$, and approximate the integral by the trapezoidal rule:

$$I[f]\approx I_{0,0}=(b-a)\frac{f(a)+f(b)}{2}=0.6051$$ (151)

We proceed with $n=1$ and $h_1=(b-a)/2=1.5$ and approximate the integral by

$$I[f]\approx I_{1,0}=\frac{h}{2}[f(a)+2f(a+h)+f(b)]=0.4968$$ (152)

We can now obtain a new method by Richardson's extrapolation to generate a more accurate result:

$$I_{1,1}=\frac{4I_{1,0}-I_{0,0}}{3}=0.46072$$ (153)

This process is carried out further recursively to generate the results shown in the table below. Every time the row index $m$ is increased by 1, the step size $h_m=h/2^m$ is halved; and every time the column index $n$ is increased by 1, a new method is obtained by more more level of recursion. All results in the nth column are generated by the same method $I_n$ but based on different step sizes; all results in the mth row are generated by different methods but all based the same step size $h_n=h/2^n$. Each element $I_{m,n}$ in the table is obtained based on its left neighbor $I_{m,n-1}$ and its top-left neighbor $I_{m-1,n-1}$ by the recursion:

$$I_{m,n}=\frac{4^nI_{m,n-1}-I_{m-1,n-1}}{4^n-1}$$ (154)

The most accurate result at each recursion level is on the diagonal of the table. In this specific case, $I_{4,4}=0.498650193$ is the most accurate result achieved after $n=4$ levels of recursion using a step size $h_4=(b-a)/2^4=(b-a)/16$.

The last row shows the relative errors of the five methods $I_{n,n}, \;n=0,\cdots,4$ along the diagonal.

$$\begin{tabular}{c\vert c\vert\vert c\vert c\vert c\vert c\vert c}... ...r&1.0641e-01&3.7928e-02&2.3465e-03&3.4978e-05&9.0721e-08\\ \hline \end{tabular}$$

The Matlab code used to generate this table is shown below:

    h=b-a;                     % initial step size
    n=5;                       % number of recursions
    m=n;
    I=zeros(n);                % initializing array I
    for i=1:m                  % for all rows
        n=2^(i-1);
        I(i,1)=(f(a)+f(b))/2;  
        for k=1:n-1
            I(i,1)=I(i,1)+f(a+k*h);
        end
        I(i,1)=h*I(i,1);       % composite trapezoidal rule
        for j=2:i              % for all columns (Richardson extrapolation)
            I(i,j)=(4^(j-1)*I(i,j-1)-I(i-1,j-1))/(4^(j-1)-1);
        end
        h=h/2;
    end