Finite Differences¶
Let $f(x) = \arctan(x)$ and $c = \sqrt{3}$, so that $f'(c) = 1/(1 + c^2) = 1/4$. We approximate $f'(c)$ by the forward, backward, and centered difference quotients $$ F(h) = \frac{f(c+h) - f(c)}{h}, \qquad B(h) = \frac{f(c) - f(c-h)}{h}, \qquad C(h) = \frac{f(c+h) - f(c-h)}{2h}, $$ for $h = 2^{-n}$, $n = 0, 1, 2, \ldots$
import numpy as np
import matplotlib.pyplot as plt
def f(x):
return np.arctan(x)
c = np.sqrt(3)
exact = 1 / 4
def forward(h):
return (f(c + h) - f(c)) / h
def backward(h):
return (f(c) - f(c - h)) / h
def centered(h):
return (f(c + h) - f(c - h)) / (2 * h)
Errors and rates¶
If $E(h) \approx C h^p$, then $E(h) / E(h/2) \approx 2^p$, so we estimate the order $p$ by the rate $$ p \approx \log_2 \frac{E(h)}{E(h/2)}. $$ We expect $p = 1$ for the forward and backward differences, and $p = 2$ for the centered difference.
methods = [forward, backward, centered]
print(f"{'n':>3} {'h':>10} | {'forward':>10} {'rate':>5} | {'backward':>10} {'rate':>5} | {'centered':>10} {'rate':>5}")
print("-" * 71)
prev = None
for n in range(11):
h = 2.0**-n
errors = [abs(exact - q(h)) for q in methods]
row = f"{n:3d} {h:10.3e}"
for i, e in enumerate(errors):
rate = f"{np.log2(prev[i] / e):5.2f}" if prev else " " * 5
row += f" | {e:10.3e} {rate}"
print(row)
prev = errors
n h | forward rate | backward rate | centered rate ----------------------------------------------------------------------- 0 1.000e+00 | 7.728e-02 | 1.653e-01 | 4.400e-02 1 5.000e-01 | 4.521e-02 0.77 | 6.642e-02 1.32 | 1.060e-02 2.05 2 2.500e-01 | 2.466e-02 0.87 | 2.989e-02 1.15 | 2.616e-03 2.02 3 1.250e-01 | 1.291e-02 0.93 | 1.421e-02 1.07 | 6.518e-04 2.00 4 6.250e-02 | 6.606e-03 0.97 | 6.932e-03 1.04 | 1.628e-04 2.00 5 3.125e-02 | 3.343e-03 0.98 | 3.424e-03 1.02 | 4.069e-05 2.00 6 1.562e-02 | 1.681e-03 0.99 | 1.702e-03 1.01 | 1.017e-05 2.00 7 7.812e-03 | 8.432e-04 1.00 | 8.483e-04 1.00 | 2.543e-06 2.00 8 3.906e-03 | 4.222e-04 1.00 | 4.235e-04 1.00 | 6.358e-07 2.00 9 1.953e-03 | 2.113e-04 1.00 | 2.116e-04 1.00 | 1.589e-07 2.00 10 9.766e-04 | 1.057e-04 1.00 | 1.058e-04 1.00 | 3.974e-08 2.00
Log-log plot¶
If $E(h) \approx C h^p$, then $\log E(h) \approx \log C + p \log h$, so on a log-log plot the errors lie on straight lines whose slope is the order $p$. The dashed and dotted lines are reference lines with slopes 1 and 2.
hs = 2.0 ** -np.arange(11)
for q, marker in zip(methods, ["o", "x", "s"]):
plt.loglog(hs, [abs(exact - q(h)) for h in hs], marker=marker, label=q.__name__)
plt.loglog(hs, 0.2 * hs, "k--", label="$O(h)$")
plt.loglog(hs, 0.02 * hs**2, "k:", label="$O(h^2)$")
plt.xlabel("$h$")
plt.ylabel("Error")
plt.legend()
plt.show()
