Floating Point Arithmetic¶
Let $x = e^h$ with $h = 2^{-n}$, and $y = 1$. We compare the computed difference $\operatorname{fl}(x) -_{\operatorname{fl}} y$ with $\operatorname{fl}(x - y)$, the exact difference rounded to the nearest double.
In [1]:
import math
Caution: math.exp and math.expm1 might not do the correct rounding, so these values may differ slightly from the lecture notes.
In [2]:
print(f"{'n':<3} | {'exp(h) - 1':<23} | {'math.expm1(h)':<23} | {'rel. error':<10}")
print("-" * 68)
for n in [*range(26, 33), None, *range(50, 55)]:
if n is None:
print(f"{':':<3} | {':':<23} | {':':<23} | {':':<10}")
continue
h = 2.0**-n
naive = math.exp(2**(-n)) - 1
accurate = math.expm1(h)
rel = abs(naive - accurate) / abs(accurate)
print(f"{n:<3d} | {naive:<23.16e} | {accurate:<23.16e} | {rel:<10.3e}")
n | exp(h) - 1 | math.expm1(h) | rel. error -------------------------------------------------------------------- 26 | 1.4901161193847656e-08 | 1.4901161304869959e-08 | 7.451e-09 27 | 7.4505805969238281e-09 | 7.4505806246794037e-09 | 3.725e-09 28 | 3.7252902984619141e-09 | 3.7252903054008080e-09 | 1.863e-09 29 | 1.8626451492309570e-09 | 1.8626451509656805e-09 | 9.313e-10 30 | 9.3132257461547852e-10 | 9.3132257504915938e-10 | 4.657e-10 31 | 4.6566128730773926e-10 | 4.6566128741615948e-10 | 2.328e-10 32 | 2.3283064365386963e-10 | 2.3283064368097468e-10 | 1.164e-10 : | : | : | : 50 | 8.8817841970012523e-16 | 8.8817841970012563e-16 | 4.441e-16 51 | 4.4408920985006262e-16 | 4.4408920985006271e-16 | 2.220e-16 52 | 2.2204460492503131e-16 | 2.2204460492503136e-16 | 2.220e-16 53 | 0.0000000000000000e+00 | 1.1102230246251565e-16 | 1.000e+00 54 | 0.0000000000000000e+00 | 5.5511151231257827e-17 | 1.000e+00
In [3]:
def forward_diff_exp(h, c=0.0):
return (math.exp(c + h) - math.exp(c)) / h
In [4]:
print(f"{'n':<3} | {'approximation':<23} | {'rel. error':<10}")
print("-" * 42)
for n in [*range(22, 28), None, *range(52, 55)]:
if n is None:
print(f"{':':<3} | {':':<23} | {':':<10}")
continue
h = 2.0**-n
approx = forward_diff_exp(h)
rel = abs(approx - 1)
print(f"{n:<3d} | {approx:<23.16e} | {rel:<10.3e}")
n | approximation | rel. error ------------------------------------------ 22 | 1.0000001192092896e+00 | 1.192e-07 23 | 1.0000000596046448e+00 | 5.960e-08 24 | 1.0000000298023224e+00 | 2.980e-08 25 | 1.0000000149011612e+00 | 1.490e-08 26 | 1.0000000000000000e+00 | 0.000e+00 27 | 1.0000000000000000e+00 | 0.000e+00 : | : | : 52 | 1.0000000000000000e+00 | 0.000e+00 53 | 0.0000000000000000e+00 | 1.000e+00 54 | 0.0000000000000000e+00 | 1.000e+00
Any differences from the lecture notes (e.g. for $n = 26$ or $n = 53$) are because math.exp is not guaranteed to round correctly.
With $c = 1$, so that $f'(1) = e$:
In [5]:
print(f"{'n':<3} | {'approximation':<23} | {'rel. error':<10}")
print("-" * 42)
for n in [*range(16, 41, 2), None, *range(50, 54)]:
if n is None:
print(f"{':':<3} | {':':<23} | {':':<10}")
continue
h = 2.0**-n
approx = forward_diff_exp(h, c=1.0)
rel = abs(approx - math.e) / math.e
print(f"{n:<3d} | {approx:<23.16e} | {rel:<10.3e}")
n | approximation | rel. error ------------------------------------------ 16 | 2.7183025674312375e+00 | 7.629e-06 18 | 2.7182870132382959e+00 | 1.907e-06 20 | 2.7182831247337162e+00 | 4.769e-07 22 | 2.7182821538299322e+00 | 1.197e-07 24 | 2.7182819098234177e+00 | 2.993e-08 26 | 2.7182818651199341e+00 | 1.349e-08 28 | 2.7182818651199341e+00 | 1.349e-08 30 | 2.7182822227478027e+00 | 1.451e-07 32 | 2.7182826995849609e+00 | 3.205e-07 34 | 2.7182846069335938e+00 | 1.022e-06 36 | 2.7182922363281250e+00 | 3.829e-06 38 | 2.7182617187500000e+00 | 7.398e-06 40 | 2.7182617187500000e+00 | 7.398e-06 : | : | : 50 | 3.0000000000000000e+00 | 1.036e-01 51 | 3.0000000000000000e+00 | 1.036e-01 52 | 4.0000000000000000e+00 | 4.715e-01 53 | 0.0000000000000000e+00 | 1.000e+00
Centered difference quotient¶
The centered difference approximation of $f'(1) = e$ is $$ f'(c) \approx \frac{f(c + h) - f(c - h)}{2h}. $$
In [6]:
def centered_diff_exp(h, c=0.0):
return (math.exp(c + h) - math.exp(c - h)) / (2 * h)
In [7]:
print(f"{'n':<3} | {'approximation':<23} | {'rel. error':<10}")
print("-" * 42)
for n in [*range(8, 35, 2), None, *range(50, 54)]:
if n is None:
print(f"{':':<3} | {':':<23} | {':':<10}")
continue
h = 2.0**-n
approx = centered_diff_exp(h, c=1.0)
rel = abs(approx - math.e) / math.e
print(f"{n:<3d} | {approx:<23.16e} | {rel:<10.3e}")
n | approximation | rel. error ------------------------------------------ 8 | 2.7182887414124934e+00 | 2.543e-06 10 | 2.7182822605182082e+00 | 1.589e-07 12 | 2.7182818554629193e+00 | 9.934e-09 14 | 2.7182818301480438e+00 | 6.213e-10 16 | 2.7182818285655230e+00 | 3.917e-11 18 | 2.7182818284491077e+00 | 3.656e-12 20 | 2.7182818283326924e+00 | 4.648e-11 22 | 2.7182818287983537e+00 | 1.248e-10 24 | 2.7182818278670311e+00 | 2.178e-10 26 | 2.7182818353176117e+00 | 2.523e-09 28 | 2.7182818055152893e+00 | 8.441e-09 30 | 2.7182819843292236e+00 | 5.734e-08 32 | 2.7182817459106445e+00 | 3.037e-08 34 | 2.7182807922363281e+00 | 3.812e-07 : | : | : 50 | 2.7500000000000000e+00 | 1.167e-02 51 | 2.5000000000000000e+00 | 8.030e-02 52 | 3.0000000000000000e+00 | 1.036e-01 53 | 0.0000000000000000e+00 | 1.000e+00
Complex step differentiation¶
The complex step approximation of $f'(1) = e$ is $$ f'(c) \approx \frac{\operatorname{Im} f(c + ih)}{h}, $$ which has error $\mathcal{O}(h^2)$ and involves no subtraction, so there is no cancellation in the numerator.
In [8]:
import cmath
def complex_step_exp(h, c=0.0):
return cmath.exp(c + 1j * h).imag / h
In [9]:
print(f"{'n':<3} | {'approximation':<23} | {'rel. error':<10}")
print("-" * 42)
for n in [*range(8, 35, 2), None, *range(200, 204)]:
if n is None:
print(f"{':':<3} | {':':<23} | {':':<10}")
continue
h = 2.0**-n
approx = complex_step_exp(h, c=1.0)
rel = abs(approx - math.e) / math.e
print(f"{n:<3d} | {approx:<23.16e} | {rel:<10.3e}")
n | approximation | rel. error ------------------------------------------ 8 | 2.7182749155161470e+00 | 2.543e-06 10 | 2.7182813963998051e+00 | 1.589e-07 12 | 2.7182818014553414e+00 | 9.934e-09 14 | 2.7182818267713138e+00 | 6.209e-10 16 | 2.7182818283535619e+00 | 3.881e-11 18 | 2.7182818284524526e+00 | 2.425e-12 20 | 2.7182818284586330e+00 | 1.516e-13 22 | 2.7182818284590193e+00 | 9.476e-15 24 | 2.7182818284590438e+00 | 4.901e-16 26 | 2.7182818284590451e+00 | 0.000e+00 28 | 2.7182818284590451e+00 | 0.000e+00 30 | 2.7182818284590451e+00 | 0.000e+00 32 | 2.7182818284590451e+00 | 0.000e+00 34 | 2.7182818284590451e+00 | 0.000e+00 : | : | : 200 | 2.7182818284590451e+00 | 0.000e+00 201 | 2.7182818284590451e+00 | 0.000e+00 202 | 2.7182818284590451e+00 | 0.000e+00 203 | 2.7182818284590451e+00 | 0.000e+00
Gaussian Elimination¶
For $a = 2^{-54}$, we solve $$ \begin{pmatrix} a & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} 1 \\ 2 \end{pmatrix} $$ by elimination and back substitution, as in the lecture: $y = (2 - 1/a) / (1 - 1/a)$ and $x = (1 - y)/a$. The exact solution is $x = 1/(1 - a)$, $y = (1 - 2a)/(1 - a)$, both very nearly 1.
In [10]:
a = 2.0**-54
x_exact = 1 / (1 - a)
y_exact = (1 - 2 * a) / (1 - a)
y = (2 - 1 / a) / (1 - 1 / a)
x = (1 - y) / a
print(f"x = {x:<23.16e} rel. error = {abs(x - x_exact) / x_exact:.3e}")
print(f"y = {y:<23.16e} rel. error = {abs(y - y_exact) / y_exact:.3e}")
x = 2.0000000000000000e+00 rel. error = 1.000e+00 y = 9.9999999999999989e-01 rel. error = 0.000e+00
np.linalg.solve uses partial pivoting (swapping the two rows, so that the very small number $a$ is not used as a pivot), which avoids the cancellation:
In [11]:
import numpy as np
A = np.array([[a, 1.0], [1.0, 1.0]])
b = np.array([1.0, 2.0])
x, y = np.linalg.solve(A, b)
print(f"x = {x:<23.16e} rel. error = {abs(x - x_exact) / x_exact:.3e}")
print(f"y = {y:<23.16e} rel. error = {abs(y - y_exact) / y_exact:.3e}")
x = 1.0000000000000000e+00 rel. error = 0.000e+00
y = 9.9999999999999989e-01 rel. error = 0.000e+00