Richardson Extrapolation¶
Let $f(x) = \arctan(x)$ and $c = 1$, so that $L = f'(c) = 1/(1 + c^2) = 1/2$. We start from the centered difference approximation $$ \phi(h) = \frac{f(c+h) - f(c-h)}{2h}, $$ whose error expansion contains only even powers of $h$. Richardson extrapolation combines values of $\phi$ at $h, h/2, h/4, \ldots$ into more accurate approximations, which we arrange in a table: $$ R(m, 0) = \phi(2^{-m} h), \qquad R(m, n) = \frac{2^{2n} R(m, n-1) - R(m-1, n-1)}{2^{2n} - 1} \quad \text{for } 1 \leq n \leq m. $$ Each column $n$ is an $O(h^{2n+2})$ approximation of $L$.
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return np.arctan(x)
c = 1.0
exact = 1 / 2
def phi(h):
return (f(c + h) - f(c - h)) / (2 * h)
The Richardson table¶
We take $h = 1$ and $N = 5$, so the first column uses $\phi(1), \phi(1/2), \ldots, \phi(1/32)$. Only the first column calls $\phi$; every other entry comes from the recursion, using the entry to its left and the entry diagonally above it.
h = 1.0
N = 5
R = np.zeros((N + 1, N + 1))
for m in range(N + 1):
R[m, 0] = phi(h / 2**m)
for n in range(1, m + 1):
R[m, n] = (2**(2 * n) * R[m, n - 1] - R[m - 1, n - 1]) / (2**(2 * n) - 1)
Errors¶
The absolute errors $E(m, n) = |L - R(m, n)|$, one row for each $m$ and one column for each $n$. The errors decrease down each column, and much faster in the later columns. The smallest errors, around $10^{-12}$, are still well above the limits of double-precision arithmetic (about $10^{-16}$), so round-off does not affect the table.
E = np.abs(exact - R)
print("m \\ n" + "".join(f"{n:>12d}" for n in range(N + 1)))
for m in range(N + 1):
print(f"{m:5d} " + "".join(f"{E[m, n]:12.3e}" for n in range(m + 1)))
m \ n 0 1 2 3 4 5
0 5.357e-02
1 1.915e-02 7.670e-03
2 5.109e-03 4.294e-04 5.335e-05
3 1.296e-03 2.508e-05 1.875e-06 1.058e-06
4 3.251e-04 1.536e-06 3.294e-08 3.706e-09 4.289e-10
5 8.136e-05 9.553e-08 5.279e-10 1.336e-11 1.116e-12 6.982e-13
Observed orders of convergence¶
Since $E(m, n) = O(h^{2(n+1)})$ and each row halves $h$, we expect $E(m-1, n) / E(m, n) \approx 2^{2(n+1)}$, so the observed order $$ p \approx \log_2 \frac{E(m-1, n)}{E(m, n)} $$ should approach $2(n+1) = 2, 4, 6, \ldots$ in columns $n = 0, 1, 2, \ldots$
Column 0 approaches 2, column 1 approaches 4, column 2 approaches 6, and column 3 is close to 8. The first entries in each column are further off, since $h$ is not yet small there. The last column has only one entry, from a fairly large $h$, so it does not show a clear pattern; increasing $N$ (more rows) would make it more apparent.
(The example in the lecture uses $c = \sqrt{2}$ instead. There, the two leading terms in the error of column 1 have opposite signs and nearly cancel for the first few values of $h$, so the observed orders in that column are irregular at first. With $c = 1$ this does not happen.)
print("m \\ n" + "".join(f"{n:>7d}" for n in range(N)))
for m in range(1, N + 1):
print(f"{m:5d} " + "".join(f"{np.log2(E[m - 1, n] / E[m, n]):7.2f}" for n in range(m)))
m \ n 0 1 2 3 4
1 1.48
2 1.91 4.16
3 1.98 4.10 4.83
4 1.99 4.03 5.83 8.16
5 2.00 4.01 5.96 8.12 8.59
Log-log plot¶
The errors of the first three columns against $h_m = 2^{-m} h$, on a log-log plot, where an error $O(h^p)$ appears as a line with slope $p$. All three are close to straight lines, with slopes about 2, 4, and 6. Column 0 goes from about $5 \times 10^{-2}$ at $h = 1$ down to about $8 \times 10^{-5}$ at $h = 1/32$; column 1 from about $8 \times 10^{-3}$ at $h = 1/2$ to $10^{-7}$; column 2 from about $5 \times 10^{-5}$ at $h = 1/4$ to $5 \times 10^{-10}$. The later columns are both smaller and steeper, so they improve much faster as $h$ decreases.
hs = h / 2.0 ** np.arange(N + 1)
for n, marker in zip(range(3), ["o", "x", "s"]):
plt.loglog(hs[n:], E[n:, n], marker=marker, label=f"column {n}, $O(h^{{{2 * n + 2}}})$")
plt.xlabel("$h$")
plt.ylabel("Error")
plt.legend()
plt.show()
