See Computational Physics (Landau, Páez, Bordeianu), Chapter 5.1 – 5.6.
These slides include material from Computational Physics. eTextBook Python 3rd Edition. Copyright © 2012 Landau, Rubin, Páez. Used under the Creative-Commons Attribution-NonCommerical-ShareAlike 3.0 Unported License.
Taking numerical derivatives is based on the elementary definition
$$ \frac{dy(t)}{dt} := \lim_{h\rightarrow 0} \frac{y(t+h) - y(t)}{h} $$Problem $h\rightarrow 0$:
number / epsilon_m
Expand $y(t)$ in Taylor series
$$ y(t+h) = y(t) + h y^{(1)}(t) + \frac{h^2}{2!} y^{(2)}(t) + \frac{h^3}{3!} y^{(3)}(t) + \cdots $$(with the derivatives $y^{(n)} \equiv \frac{d^{n}y}{dt^n}$)
Neglect $\mathcal{O}(h^2)$ terms:
$$ \frac{dy(t)}{dt} \simeq D_\text{fd} y(t) \equiv \frac{y(t+h) - y(t)}{h} $$Image from Computational Physics. eTextBook Python 3rd Edition. Copyright © 2012 Landau, Rubin, Páez. Used under the Creative-Commons Attribution-NonCommerical-ShareAlike 3.0 Unported License.
Substitute the Taylor series into $D_\text{fd} y(t)$
$$ D_\text{fd} y(t) = y^{(1)}(t) - \frac{h}{2} y^{(2)}(t) + \cdots $$so the error is first order in $h$: $$ y^{(1)}(t) - D_\text{fd} y(t) = \frac{1}{2} h y^{(2)}(t) + \mathcal{O}(h^2) $$
Quadratic function $$ y(t) = a + bt^2 \quad \text{and}\quad y^{(1)}(t) = 2bt $$ Forward difference approximation $$ D_\text{fd} y(t) = \frac{y(t+h) - y(t)}{h} = 2bt + bh $$ is only good for small $h \ll \frac{1}{b}$.
Improve over the (rather bad) error in the forward difference algorithm by using the central difference.
Image from Computational Physics. eTextBook Python 3rd Edition. Copyright © 2012 Landau, Rubin, Páez. Used under the Creative-Commons Attribution-NonCommerical-ShareAlike 3.0 Unported License.
Applied to the parabola $y(t) = a + bt^2$ yields the exact result: $$ D_\text{cd} = 2bt $$
Insert Taylor expansion...
$$ y^{(1)}(t) - D_\text{cd} y(t) = \frac{1}{24} h^2 y^{(3)}(t) + \mathcal{O}(h^4) $$The error is second order in $h$.
This is generally much better than the forward difference.
"Generally" means "when $y(t)$ is sufficiently smooth", i.e. \begin{align} y^{(3)} \frac{h^2}{24} &\ll y^{(2)} \frac{h}{2} \quad\text{or}\\ 12\frac{y^{(2)}}{y^{(3)}} &\gg h. \end{align}
Other difference algorithms exist with even better error behavior (see homework) but all have problems with data that are not smooth. In this case you might have to smooth the data in advance by interpolation or kernel density estimates (KDE).
Round-off error of differentiation: difference of two similar numbers cannot be better than machine precision, always divided by small number $h$:
$$ \epsilon_\text{ro} \approx \frac{\epsilon_m}{h} $$Forward difference algorithmic error $$ \epsilon_\text{app}^\text{fd} = y^{(2)} \frac{h}{2} $$
What is the most accurate value of $h$ that we can choose for $D_\text{fd}$?
Cross-over from approximation error to round-off error around
$$ \epsilon_\text{app} = \epsilon_\text{ro} $$For the forward difference $$ y^{(2)} \frac{h_\text{fd}}{2} = \frac{\epsilon_m}{h_\text{fd}}\\ h_\text{fd}^2 = \frac{2\epsilon_m}{y^{(2)}} $$
For rough estimate, assume $y^{(1)} \approx y^{(2)} \approx y^{(3)} \approx 1$:
import numpy as np
def h_fd(eps_m, y2=1):
return np.sqrt(2*eps_m/y2)
def eps_fd(eps_m, y2=1.):
h = h_fd(eps_m, y2=y2)
return 2*eps_m/h
eps_m = 1e-15
print("h_fd = {0:.2e}".format(h_fd(eps_m, y2=1)))
print("eps_fd = {0:.2e}".format(eps_fd(eps_m, y2=1)))
h_fd = 4.47e-08 eps_fd = 4.47e-08
For double precision ($\epsilon_m = 10^{-15}$) $$ h_\text{fd} \approx 4\times 10^{-8}\\ \epsilon_\text{tot}^\text{fd} \approx 2\epsilon_\text{ro} = \frac{2\epsilon_m}{h_\text{fd}} = h_\text{fd} = 4 \times 10^{-8} $$
(Note: simplified because $y^{(2)} \approx 1$.)
Do you expect the "best" $h_\text{cd}$ for the central difference algorithm to be larger or smaller than $h_\text{fd}$?
Generally $h_\text{cd} > h_\text{fd}$ because the better algorithm allows you to use a bigger step size for a smaller error.
and also the extended difference algorithm
\begin{align} D_\text{ep} y(t) &\equiv \frac{4 D_\text{cd}y(t, h/2) - D_\text{cd}y(t, h)}{3} \\ &= \frac{8\big(y(t+h/4) - y(t-h/4)\big) - \big(y(t+h/2) - y(t-h/2)\big)}{3h} \end{align}def D_fd(y, t, h):
"""Forward difference"""
return (y(t + h) - y(t))/h
def D_cd(y, t, h):
"""Central difference"""
# implement
def D_ed(y, t, h):
"""Extended difference"""
# implement
Test function: $y(t) = \cos t$
Try to do the above for all three algorithms
git pull
the course resources and work with the notebook 06_differentiation/06_differentiation-students.ipynb.
def D_fd(y, t, h):
"""Forward difference"""
return (y(t + h) - y(t))/h
def D_cd(y, t, h):
"""Central difference"""
return (y(t + h/2) - y(t - h/2))/h
def D_ed(y, t, h):
"""Extended difference"""
return (8*(y(t + h/4) - y(t - h/4)) - (y(t + h/2) - y(t - h/2)))/(3*h)
import numpy as np
# test function: np.cos
# Analytical derivative
def y2(t):
return -np.sin(t)
t_values = np.array([0.1, 1, 100], dtype=np.float64)
Use numpy functions for everything because then you can operate on all t_values
at once.
Note that we pass a function y
to the forward difference function D_fd
and we can also pass a whole array of t_values
!
D_fd(np.cos, t_values, 0.1)
array([-0.14937587, -0.86706184, 0.46244209])
D_fd(np.cos, t_values, 1e-12)
array([-0.09980905, -0.84154905, 0.50359716])
Compute the exact derivatives (again, operate on all $t$ together... start thinking in numpy arrays!)
y2(t_values)
array([-0.09983342, -0.84147098, 0.50636564])
Calculation of the absolute error: subtract the two arrays that you got previously:
D_fd(np.cos, t_values, 0.1) - y2(t_values)
array([-0.04954246, -0.02559086, -0.04392355])
def error(Dxx, y, y2, t, h):
"""Relative error.
Note: will fail for exact 0.
"""
y2_val = y2(t)
# return (Dxx(y, t, h) - y2_val/y2_val
return Dxx(y, t, h)/y2_val - 1
Note that we pass again a general function for the difference operator so that we can use error()
with D_fd()
, D_cd()
and D_ep()
.
error(D_fd, np.cos, y2, t_values, 0.1)
array([ 0.49625125, 0.03041205, -0.08674275])
error(D_fd, np.cos, y2, t_values, 1e-12)
array([ -2.44073917e-04, 9.27754603e-05, -5.46734793e-03])
Plot $\log_{10} |E(h)|$ against $\log_{10} h$.
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
#matplotlib.style.use('ggplot') # check out available styles with matplotlib.style.available
matplotlib.style.use('seaborn-talk')
h_values = 10**(np.arange(-15, -1, 0.1))
abs_errors = np.abs(error(D_fd, np.cos, y2, 0.1, h_values))
plt.loglog(h_values, abs_errors, label="t=0.1")
[<matplotlib.lines.Line2D at 0x11cce0da0>]
Plot the three different $t$ values together in one plot:
for t in t_values:
abs_errors = np.abs(error(D_fd, np.cos, y2, t, h_values))
plt.loglog(h_values, abs_errors, label=r"$t={}$".format(t))
ax = plt.gca()
ax.legend(loc="best")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$|E_\mathrm{fd}|$")
<matplotlib.text.Text at 0x11cbdcac8>
Plot the three different $t$ values together in one plot:
for t in t_values:
abs_errors = np.abs(error(D_cd, np.cos, y2, t, h_values))
plt.loglog(h_values, abs_errors, label=r"$t={}$".format(t))
ax = plt.gca()
ax.legend(loc="best")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$|E_\mathrm{cd}|$")
<matplotlib.text.Text at 0x11cec9cc0>
Plot the three different $t$ values together in one plot:
for t in t_values:
abs_errors = np.abs(error(D_ed, np.cos, y2, t, h_values))
plt.loglog(h_values, abs_errors, label=r"$t={}$".format(t))
ax = plt.gca()
ax.legend(loc="best")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$|E_\mathrm{ed}|$")
<matplotlib.text.Text at 0x11cfbebe0>