See Computational Physics (Landau, Páez, Bordeianu), Chapter 5.5
Forward difference and central difference (discussed in class):
\begin{align} D_\text{fd} y(t) &\equiv \frac{y(t+h) - y(t)}{h} \\ D_\text{cd} y(t) &\equiv \frac{y\Big(t + \frac{h}{2}\Big) - y\Big(t - \frac{h}{2}\Big)}{h} \end{align}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_ep(y, t, h):
"""Extended difference"""
# implement
Test function: $y(t) = \cos t$
Try to do the above for all three algorithms
import numpy as np
# test function
y = np.cos
# Analytical derivative (use np.cos, np.sin, np.exp or whatever else you need)
def y2(t):
# IMPLEMENT ME (remove `pass`)
pass
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(y, t_values, 0.1)
Compute the exact derivatives (again, operate on all $t$ together... start thinking in numpy arrays!)
y2(t_values)
Calculation of the absolute error: subtract the two arrays that you got previously:
# IMPLEMENT ME
def error(Dxx, y, y2, t, h):
"""Relative error."""
# INCOMPLETE ... replace None with appropriate terms
return (Dxx(y, t, h) - None)/None
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, y, y2, t_values, 0.1)
Plot $\log_{10} |E(h)|$ against $\log_{10} h$.
import matplotlib.pyplot as plt
%matplotlib inline
h_values = 10**(np.arange(-15, -1, 0.1))
abs_errors = np.abs(error(D_fd, y, y2, 0.1, h_values))
# INCOMPLETE, replace None with appropriate terms
plt.loglog(None, None, label="t=0.1")
Plot the three different $t$ values together in one plot:
# INCOMPLETE: replace None values by appropriate terms
for t in t_values:
abs_errors = None
plt.loglog(None, None, label=r"$t={}$".format(t))
ax = plt.gca()
ax.legend(loc="best")
ax.set_xlabel(r"$h$")
ax.set_ylabel(r"$|E|$")