Differentiation by Richardson Extrapolation

The Richardson extrapolation can be applied to improve numerical differentiation. Recall that the derivative $D[f]=f'(x)$ can be approximated by the central difference based on two neighboring points $x-h$ and $x+h$ with step size $h$ and a truncation error term $O(h^2)$, as shown in Eq. (6):

$$D[f] = \frac{f(x+h)-f(x-h)}{2h} -\left[\frac{f'''(x)}{3!}h^2+\frac{f^{(5)}(x)}{5!}h^4+\cdots\right]$$

$$= D_0(h)+C_1h^2+C_2h^4+C_3h^6+C_4h^8+\cdots=D_0(h)+C_2h^2+O(h^4)$$ (125)

where

$$D_0(h)=\frac{f(x+h)-f(x-h)}{2h}$$ (126)

The accuracy of this method $D_0(h)$ can be improved by the Richardson extrapolation considered previously. We first replace $h$ in the equation above by $h/2$:

$$D[f]=D_0(h/2)+C_2h^2/4+O(h^4)$$ (127)

and then we get

$$D[f]=\frac{4D_0(h/2)-D_0(h)}{3}+O(h^4)=D_1(h/2)+O(h^4)$$ (128)

where

$$D_1(h/2) = \frac{4D_0(h/2)-D_0(h)}{3} =\frac{1}{3}\left[ 4\frac{f(x+h/2)-f(x-h/2)}{h} -\frac{f(x+h)-f(x-h)}{2h}\right]$$

$$= \frac{1}{6h}\left[f(x-h)-8f(x-h/2)+8f(x+h/2)-f(x+h)\right]$$ (129)

is a new method based on the four points $x-h,\;x-h/2,\;x+h/2$ and $x+h$ 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 $D_1(h/2)$ by $h/2$ and repeating the same steps we get:

$$D[f] = \frac{16 D_1(h/4)-D_1(h/2)}{15}+O(h^6)=D_2(h/4)+O(h^6)$$ (130)

where we defined a new method:

$$D_2(h/4) = \frac{16 D_1(h/4)-D_1(h/2)}{15}$$

$$= \frac{1}{15}\;\frac{16}{3h}[f(x-h/2)-8f(x-h/4) +8f(x+h/4)+f(x+h/2)]$$

$$-\frac{1}{15}\;\frac{1}{6h}[f(x-h)-8f(x-h/2)+8f(x+h/2)-f(x+h)]$$

$$= \frac{1}{90h}[-f(x-h)+40f(x-h/2)-256f(x-h/4)$$

$$+256f(x+h/4)-40f(x+h/2)+f(x+h)]$$ (131)

based on six points $x\pm h,\;x\pm h/2,\;x\pm h/4$, with a step size $h/4$ and an error term $O(h^6)$. Note that $D_1(h)$ and $D_2(h)$ are actually the same methods previously derived by the undetermined coefficients approach. However, we now realize that by the process discussed above can be recursively carried out to obtain methods $D_3$, $D_4$, etc., with progressively more accurate results.

The process of the Richardson extrapolation can be repeatedly carried out to improve the accuracy.

Example:

Find the derivative of the Gaussian function $f(x)=e^{-x^2}$ at $x=1$. Initially, we set $n=m=0$ and use the divided difference method with a step size $h_0=h=1$:

$$f'(x)\approx D_{0,0}=\frac{f(x+h)-f(x-h)}{2h}=-0.4908$$ (132)

We proceed with $n=0$ for the same method but $m=1$ for a smaller step side $h_1=h_0/2^m=h_0/2=0.5$:

$$D_{1,0}=\frac{f(x+h/2)-f(x-h/2)}{2}=-0.6734$$ (133)

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

$$D_{1,1}=\frac{4D_{1,0}-D_{0,0}}{3}=-0.73425$$ (134)

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 one more level of recursion. All results in the nth column are generated by the same method $D_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 $D_{m,n}$ in the table is obtained based on its left neighbor $D_{m,n-1}$ and its top-left neighbor $D_{m-1,n-1}$ by the recursion:

$$D_{m,n}=\frac{4^nD_{m,n-1}-D_{m-1,n-1}}{4^n-1}$$ (135)

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

The last row shows the relative errors of the five methods $D_{n,n}, \;n=0,\cdots,4$ along the diagonal compared with the ground truth

$$f'(1)=\frac{d}{dx}e^{-x^2}\bigg\vert _{x=1}=-2x\, e^{-x^2}\bigg\vert _{x=1}=-0.73575888$$ (136)

$$\begin{tabular}{c\vert c\vert\vert c\vert c\vert c\vert c\vert c}... ...r&2.4492e-01&1.5042e-03&3.4678e-04&2.0553e-06&1.6927e-09\\ \hline \end{tabular}$$

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

    h=1;                                    % inital step size
    n=5;                                    % number of recursions
    m=n;
    D=zeros(n);                             % initializing array D
        for i=1:m                           % for all rows 
            n=2^(i-1);
            D(i,1)=(f(x+h)-f(x-h))/2/h;
            for j=2:i                       % for all columns
                D(i,j)=(4^(j-1)*D(i,j-1)-D(i-1,j-1))/(4^(j-1)-1);
            end
            h=h/2;
        end