Extrapolation methods

Extrapolation methods exploit the error estimate to compute better approximations. Such methods can be used for any approximation methods but here we illustrate it for numerical integration.

1Aitken extrapolation¶

Suppose we have an asymptotic error formula

I - I_n \approx \frac{c}{n^p}, \qquad p > 0

(1)

If we know $I$ then we can estimate $p$ from

\frac{I - I_n}{I - I_{2n}} = 2^p

(2)

How can we estimate the convergence rate $p$ when we dont know $I$ ? First compute $I_n$ , $I_{2n}$ , $I_{4n}$ . Then

\begin{aligned} \frac{I_{2n} - I_n}{I_{4n} - I_{2n}} &= \frac{ (I-I_n) - (I - I_{2n})}{(I - I_{2n}) - (I - I_{4n})} \\ &= \frac{ c/n^p - c/(2n)^p}{c/(2n)^p - c/(4n)^p} \\ &= 2^p \end{aligned}

(3)

Hence, the convergence rate is

p = \frac{ \log\frac{I_{2n} - I_n}{I_{4n} - I_{2n}}}{\log 2}

(4)

We can obtain a more accurate estimate of $I$ as follows. Again using the error estimate

\frac{I - I_n}{I - I_{2n}} = 2^p = \frac{I - I_{2n}}{I - I_{4n}}

(5)

(I - I_n)(I - I_{4n}) = (I - I_{2n})^2

(6)

Solving this

I \approx \tilde I_{4n} := I_{4n} - \frac{(I_{4n} - I_{2n})^2}{(I_{4n} - I_{2n}) - (I_{2n} - I_n)}

(7)

Then $\tilde I_{4n}$ is more accurate than $I_{4n}$ .

1.1Example¶

Approximate

I = \int_0^1 x \sqrt{x} \ud x

(8)

using Simpson rule with $n=16,32,64$ intervals (we have $n+1$ points). Then

I - I_{64} \approx -4.29 \times 10^{-7}, \qquad I - \tilde I_{64} \approx 6.13 \times 10^{-10}

(9)

Clearly $\tilde I_{64}$ is more accurate than $I_{64}$ . Moreover

\tilde I_{64} - I_{64} \approx 0.399999999387 - 0.400000429413 = -4.30026 \times 10^{-7} \approx I - I_{64}

(10)

1.2Remark¶

(a) If we have computed $I_{4n}$ then we can compute $I_{2n}$ and $I_n$ with very little effort since the same function values are used^[1]. (b) We can then also obtain a more accurate estimate $\tilde I_{4n}$ very easily. (c) We also have an error estimate

\tilde I_{4n} - I_{4n} \approx I - I_{4n}

(11)

If $|\tilde I_{4n} - I_{4n}| \le \epsilon$ , then we can assume that $|I - \tilde I_{4n}| \le |I - I_{4n}| \le \epsilon$ and use $\tilde I_{4n}$ as the estimate of the integral.

2Richardson extrapolation¶

This is a more sophisticated way of performing the extrapolation. Let us illustrate this technique using the trapezoidal method for which we know the error estimate of the form

I - I_n = \frac{d_2^{(0)}}{n^2} + \frac{d_4^{(0)}}{n^4} + \ldots + \frac{d_{2m}^{(0)}}{n^{2m}} + F_{n,m}

(12)

when $f \in \cts^{2m+2}[a,b]$ . For $n$ even

I - I_{n/2} = \frac{4 d_2^{(0)}}{n^2} + \frac{16 d_4^{(0)}}{n^4} + \frac{64 d_6^{(0)}}{n^6} + \ldots

(13)

Then

4(I - I_n) - (I - I_{n/2}) = -\frac{12 d_4^{(0)}}{n^4} - \frac{60 d_6^{(0)}}{n^6} + \ldots

(14)

so that

I = \frac{1}{3}(4 I_n - I_{n/2}) -\frac{4 d_4^{(0)}}{n^4} - \frac{20 d_6^{(0)}}{n^6} + \ldots

(15)

Define

I_n^{(1)} = \frac{4}{3} I_n^{(0)} - \frac{1}{3} I_{n/2}^{(0)}, \qquad n \ge 2 \textrm{ and $n$ is even}

(16)

with

I_n^{(0)} = I_n

(17)

$I_n^{(1)}$ is called the Richardson extrapolation of $I_n$ . Note that $I_n^{(1)}$ uses the same data as $I_n$ but has a convergence rate of 4 while $I_n$ converges at rate 2. But what is $I_n^{(1)}$ ?

\begin{aligned} I_n^{(1)} &= \frac{4}{3} h [\shalf f_0 + f_1 + f_2 + \ldots + f_{n-1} + \shalf f_n] \\ & \quad - \frac{1}{3} (2h)[\shalf f_0 + f_2 + f_4 + \ldots + f_{n-2} + \shalf f_n] \\ &= \frac{h}{3}[f_0 + 4 f_1 + 2 f_2 + \ldots + 4 f_{n-2} + 2 f_{n-2} + f_n ] \end{aligned}

(18)

which is just Simpson rule !!!

This procedure can be continued further. Let $n$ be a multiple of 4. We have already shown that

I - I_n^{(1)} = \frac{d_4^{(1)}}{n^4} + \frac{d_6^{(1)}}{n^6} + \ldots

(19)

where

d_1^{(1)} = -4 d_4^{(0)}, \qquad d_6^{(1)} = -20 d_6^{(0)}, \ldots

(20)

and so

I - I_{n/2}^{(1)} = \frac{16 d_4^{(1)}}{n^4} + \frac{64 d_6^{(1)}}{n^6} + \ldots

(21)

Hence

16(I - I_n^{(1)}) - (I - I_{n/2}^{(1)}) = -\frac{48 d_6^{(1)}}{n^6} + \ldots

(22)

and

I = \frac{16 I_n^{(1)} - I_{n/2}^{(1)}}{15} - \frac{48 d_6^{(1)}}{15 n^6} + \ldots

(23)

Then

I_n^{(2)} = \frac{16 I_n^{(1)} - I_{n/2}^{(1)}}{15}, \qquad n \ge 4, \quad \textrm{$n$ is multiple of 4}

(24)

converges at the rate 6 and is same as Boole’s rule.

2.1Richardson error estimate for $I_n^{(1)}$ (Simpson rule)¶

From

I = \frac{16 I_n^{(1)} - I_{n/2}^{(1)}}{15} - \frac{48 d_6^{(1)}}{15 n^6} + \ldots

(25)

we get

\begin{aligned} I - I_n^{(1)} &= \frac{16 I_n^{(1)} - I_{n/2}^{(1)}}{15} - I_n^{(1)} - \frac{48 d_6^{(1)}}{15 n^6} + \ldots \\ &= \frac{1}{15}[ I_n^{(1)} - I_{n/2}^{(1)} ] + \order{\frac{1}{n^6}} \end{aligned}

(26)

Since $I_n^{(1)} - I_{n/2}^{(1)} = \order{1/n^4}$ , we can use the error estimate

I - I_n^{(1)} \approx \frac{1}{15}[ I_n^{(1)} - I_{n/2}^{(1)} ]

(27)

The right hand side can be easily computed.

3Romberg integration¶

Using the error formula of trapezoidal rule, we can generalize Richardson extrapolation to arbitrary orders. For $n =$ some integer multiple of $2^k$

I_n^{(k)} = \frac{4^k I_n^{(k-1)} - I_{n/2}^{(k-1)}}{4^k - 1}

(28)

and

I - I_n^{(k)} = \order{\frac{1}{n^{2k+2}}}

(29)

Romberg integration is defined by

J_k(f) = I_{2^k}^{(k)}, \qquad k=0,1,2,\ldots

(30)

For example, if $k=3$ , compute

\begin{array}{cccc} I_1^{(0)} & & & \\ I_2^{(0)} & I_2^{(1)} & & \\ I_4^{(0)} & I_4^{(1)} & I_4^{(2)} & \\ I_8^{(0)} & I_8^{(1)} & I_8^{(2)} & I_8^{(3)}\\ \end{array}

(31)

The first column is trapezoid rule. The second column is generated from the first column using ([eq:romiter]{reference-type=“ref” reference=“eq:romiter”}) with $k=1$ . The third column is generated from second column and the fourth from the third. Then $I_8^{(3)}=J_3(f)$ is the final answer.

Footnotes¶

This is true for methods like trapezoid, mid-point, Simpson but is not true for Gauss rules.
↩

Numerical Analysis

Special techniques