import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.dpi'] = 120
''' Calculating derivatives in Python '''
import scipy.misc as msc
def f(x): return (np.exp(x)-1.234576)**2-1
msc.derivative(np.exp, 1)
''' Solving equations '''
import scipy.optimize as opt
opt.fsolve(f, 1)
''' ... and finding minimum '''
opt.minimize(f, 1)
Derivatives:
Our goal is to learn how to solve numerically partial differential equation (PDE) such as the Black–Scholes equation
$$\frac{\partial V}{\partial t} + \frac{1}{2}\sigma^2 S^2 \frac{\partial^2 V}{\partial S^2} + rS\frac{\partial V}{\partial S} - rV = 0.$$But before we solve anything we have to understand how the derivatives are calculated by computers.
Let us recall a definition
$$f'(x)=\lim_{h\to0}\frac{f(x+h)-f(x)}{h}=\lim_{h\to0}\frac{f(x)-f(x-h)}{h}.$$Can we use this definition to calculate the derivative numerically?
But how to choose $h$?
If we take $h$ too large we will get an inacurate value (Taylor series)
$$f(x+h)\approx f(x)+f'(x)h+\frac{1}{2!}f''(x)h^2+\ldots=\sum_{n=0}^\infty\frac{f^{(n)}(x)}{n!}h^n$$$$D_1(h)=\frac{f(x+h)-f(x)}{h}=\color{green}{f'(x)}+\color{red}{\frac{1}{2!}f''(x)h}+\frac{1}{3}f'''(x)h^2+\ldots$$Taking the average $\frac{1}{2}(D_1(h)+D_1(-h))$ we find a central difference scheme
$$D_2(h)=\frac{f(x+h)-f(x-h)}{2h}=\color{green}{f'(x)}+\color{red}{\frac{1}{6}f'''(x)h^2}+\mathcal{O}(h^4)$$Both methods require calculation of $f$ twice, but the second gives smaller error ($h\ll 1$).
Forward difference:
$$D_{2f}(h)=\frac{-3f(x)+4f(x+h)-f(x+2h)}{2h}=\color{green}{f'(x)}+\color{red}{\frac{1}{3}f'''(x)h^2}+\mathcal{O}(h^3)$$Note that central difference has twice smaller error than the forward difference, and in case of vanishing $f'''$ even the order is different.
By taking a combination of more points, we can cancel higher terms of the Taylor expansion, for example from five points we obtain a fourth order five point-stencil method.
$$D_\color{red}{4}(h)=\frac{f(x-2h)-8f(x-h)+8f(x+h)-f(x+2h)}{12h}+\mathcal{O}(h^\color{red}{4})$$Note: because of the Runge phenomenon higher order schemes for evenly separated points are unadvised.
The list of the coefficients for central and forward can be found quite easily
Second and higher derivatives can also be found in the similar way:
$$\small f''(x)\approx D^2_4(h)=\frac{-f(x-2h)+16f(x-h)-30f(x)+16f(x+h)-f(x+2h)}{12h^2}+\mathcal{O}(h^4)$$So, taking smaller $h$ always gives better results
The numbers in computer have finite precision $x\in(a-\epsilon |a|, a+\epsilon|a|)$ have the same floating point value.
$\epsilon=10^{-8}$ for floats and $10^{-16}$ for doubles.
Subtracting similar values gives huge loss of the relative accuracy
$$1.2345678(\pm10^{-8})-1.2345666(\pm10^{-8})=0.0000012(\pm10^{-8})$$but we have lost 6 significant digits!
We have two types of errors:
The best choice of $h$ minimizes the total error $\Delta=\Delta_m+\Delta_r=ah^n+bh^{-1}=\min$.
This gives rough estimation for optimal $h\sim\epsilon^{1/(n+1)}$ and $\Delta\sim\epsilon^{n/(n+1)}$
def showTheoreticalError(n):
with plt.xkcd():
fig = plt.figure(figsize=(9,4.5))
ax = fig.add_axes((0.1, 0.2, 0.8, 0.7))
ax.spines['right'].set_color('none')
ax.spines['top'].set_color('none')
title = 'Schematic error of the FDM for $\mathcal{O}$';
plt.title('{:s}$(h^{:d})$'.format(title, n) )
plt.xlabel('spaceing $h$')
plt.ylabel("error $\Delta$")
plt.xticks([])
plt.yticks([])
h = np.linspace(0.0001,5, 200)
plt.ylim(0,5)
plt.xlim(0,5)
plt.plot(h, h**n, dashes=[4,2])
plt.plot(h, 1/h, dashes=[2,2])
plt.plot(h, 1/h+h**n)
plt.legend('best')
plt.legend(['method error $\Delta_m(h)$', 'machine error $\Delta_r$', 'total error $\Delta_m+\Delta_r$'])
plt.show()
showTheoreticalError(n=1)
def D1(x, h, f): # First order FD
return (f(x+h)-f(x))/h;
def D2(x, h, f): # 2nd order FD
return (f(x+h)-f(x-h))/(2*h);
def D4(x, h, f): # 4th order FD
return (f(x-2*h)-8*f(x-h)+8*f(x+h)-f(x+2*h))/(12*h)
def f (x): # Example function
return np.sin(x);
h = np.logspace(-16, 0, num=500, endpoint=True) # scan through different h
x = 1.0;
D1f = D1(x, h, f);
df = np.cos(x); # exact f'
Error1 = np.abs(D1f-df);
Error2 = np.abs(D2(x, h, f)-df);
Error4 = np.abs(D4(x, h, f)-df);
import matplotlib.style
import matplotlib as mpl
mpl.style.use('default')
def plot_FD_errors():
plt.figure(figsize=(6,3.5), dpi=150)
plt.title('Error of the finite difference methods for $f(x)=\sin(x)$')
plt.xlabel('$h$')
plt.ylabel("$|D_nf - f'|$")
plt.grid(True)
plt.loglog(h, Error1, basex=10)
plt.loglog(h, Error2, basex=10)
plt.loglog(h, Error4, basex=10)
plt.legend('best')
plt.legend(['$\mathcal{O}(h)$', '$\mathcal{O}(h^2)$', '$\mathcal{O}(h^4)$'])
plt.show()
plot_FD_errors()
Note that in the first order scheme we lose up to 8 significant digits for the best choice of $h$! The relative error grew $10^8$ times.
Other schemes are better with the best accuracy $\Delta\sim \epsilon^{n/(n+1)}$
Can we do better: For the $\mathcal{O}(h^2)$ symmetric scheme we have
$$D_2(h) = f'+\alpha_2 h^2f'''+\alpha_4 h^4 f^{(v)} + \alpha_6 h^6f^{(vii)} + \cdots$$$$D_2(h/2) = f'+\frac{1}{4}\alpha_2 h^2f'''+\frac{1}{16}\alpha_4 h^4f^{(v)} + \frac{1}{64}\alpha_6 h^6f^{(vii)}+ \cdots$$$$\frac{4D_2(h/2)-D_2(h)}{3} = {\color{green}{f'}} - \color{red}{\frac{\alpha_4f^{(v)}h^4}{4}} +\frac{5\alpha_6 f^{(vii)}h^6}{16}+\cdots$$Using twice $\mathcal{O}(h^2)$ scheme we can obtain $\mathcal{O}(h^4)$ result.
Calculating $D_2(h/4)$ we can cancel higher terms
Let us introduce a grid: $x_n=nh, n=0,\ldots,N-1$.
Function $f(x)$ can be stored as a vector $y_n=f(x_n)$.
A vector of the first derivatives $y'_n\approx f'(x_n)$ can be calculated (using central $\mathcal{O}(h^2)$ scheme):
$$y'_n=\frac{-y_{n-1}+y_{n+1}}{2h},\qquad n>0 \wedge n<N-1.$$Note, that for $n=0$ and $n=N-1$ we would need to know $y_{-1}$ and $y_N$.
They can be found either from boundary conditions, or from symmetries (pairity, periodicity etc).
Alternatively we could use appropriate forward or backward schemes:
$$y'_0=\frac{-3y_0 + 4y_1 -y_2}{2h}$$$$y'_{N-1}=-\frac{-3y_{N-1} + 4y_{N-2} -y_{N-3}}{2h}$$All these expressions can be written in a matrix form $y'=Dy$ where $y,y'\in\mathbb{R}^{N}, D\in\mathbb{R}^{N\times N}$
$${\scriptsize D=\frac{1}{2h}\left(\begin{array}{rrrrrrrrr} -3& 4&-1\\ -1& 0& 1\\ & -1& 0& 1\\ && -1& 0& 1\\ &&&\ddots&\ddots&\ddots\\ &&&&&-1& 0& 1\\ &&&&& 1& -4& 3 \end{array} \right)}$$Remarks
Let us compare times needed to calculate $y''_n\approx\frac{1}{12h^2}(-y_{n-2}+16y_{n-1}-30y_n+16y_{n+1}-y_{n+2})$
for a whole vector without boundary points.
import time
def measure_time_Derivative2_loop(y, h, N, D):
fac = (12*h*h); start = time.time()
dy = []
for n in range(2, N-3):
dy.append( (-y[n-2]+16*y[n-1]-30*y[n]+16*y[n+1]-y[n+2])/fac );
return time.time()-start, dy
def measure_time_Derivative2_vectorize(y, h, N, D):
fac = (12*h*h); start = time.time()
dy = (-30/fac)*y[2:N-3] + (16/fac)*(y[1:N-4]+y[3:N-2]) - (1/fac)*(y[0:N-5]+y[4:N-1]);
return time.time()-start, dy
def measure_time_Derivative2_matrix(y, h, N, D):
start = time.time()
dy = D*y
return time.time()-start, dy
import scipy.sparse as sp # povides sparse matrices
def init_vectors(N):
x = np.linspace(-10, 10, N); # grid points
y = np.exp(-x**2); # example function
h = x[1]-x[0] # grid spaceing
u = np.ones(N)/(12*h*h)
data = np.array([-u, 16*u, -30*u, 16*u, -u]) # array Nx5 to fill the matrix
diags = np.array([-2, -1, 0, 1, 2]) # which diagonals to fill
D_sparse = sp.spdiags(data, diags, N, N) # sparse 5-diagonal matrix
return h, x, y, D_sparse
N=10000
h, x, y, D = init_vectors(N)
T1, dy = measure_time_Derivative2_loop(y, h, N, D)
T2, dy = measure_time_Derivative2_vectorize(y, h, N, D)
T3, dy = measure_time_Derivative2_matrix(y, h, N, D)
D_full = D.toarray(); # Full matrix from sparse D
T4, dy = measure_time_Derivative2_matrix(y, h, N, D_full)
print ("{:e} {:e} {:e} {:e}".format(T1, T2, T3, T4))
print ("Vectorization is {:.2f} times faster than loops".format(T1/T2))
print ("Vectorization is {:.2f} times faster than sparse matrix multiplication".format(T3/T2))
print ("Full matrix multiplication is {:.2f} slower than sparse".format(T4/T3))
def measure_all_times(N, iterations):
N = int(N)
h, x, y, D = init_vectors(N)
A = [N]
for f in [measure_time_Derivative2_loop, measure_time_Derivative2_vectorize, measure_time_Derivative2_matrix]:
T = []
for k in range(iterations):
t, dy = f(y, h, N, D)
T.append(t*1e6)
A.append(np.sum(np.sort(T)[2:-2])/(iterations-4))
if (N<10000): ## Full matrix is slow, so we measure it only for small N
f = measure_time_Derivative2_matrix
D_full = D.toarray()
T = []
for k in range(iterations):
t, dy = f(y, h, N, D_full)
T.append(t*1e6)
A.append(np.sum(np.sort(T)[2:-2])/(iterations-4))
else:
A.append(np.nan);
return A;
def measure_times_through_N():
measured_times = []
measure_points_number=50
for N in np.linspace(1, 5, measure_points_number):
A = measure_all_times(10**N, 15)
measured_times.append(A)
measured_times = np.reshape(measured_times, (measure_points_number,5))
np.savetxt('speed_measurements.dat', measured_times)
def show_measured_times():
import warnings
warnings.filterwarnings('ignore')
measured_times = np.loadtxt('speed_measurements.dat')
plt.figure(figsize=(9,4.5))
plt.title('Calculation time')
plt.xlabel('$N$')
plt.ylabel("Time $t [\mu s]$")
plt.grid(True)
plt.loglog(measured_times[:,0], measured_times[:,1:5], basex=10)
C = np.loadtxt('Ctests/CSpeedKing.dat', usecols=range(2))
plt.loglog(C[:,0], C[:,1], basex=10, dashes=[2,2], color='black')
plt.legend(['loop', 'vectorisation', 'sparse matrix', 'full matrix', 'C simple loop'])
plt.show()
$^*$ append reallocates list doubling its length, but preallocation is not always faster.
Standard loops have a large coefficient, some $10-100$ larger than vectorization.
# measure_times_through_N()
show_measured_times()
Note that loops in C/C++ can be faster than vectorization in Python, there is no need to create temporary array objects.
Low level language can scan memory more efficiently and can use registers (eg. BLAS library).
The function $\equiv$ values at certain points $y_i=f(x_i)$ but we can choose different parametrization (basis):
Polynomial base $\displaystyle f(x)=\sum_{n=0}^Na_nx^n$ $$f'(x)=\sum_{n=0}^Nn a_nx^{n-1}=\sum_{n=0}^{N-1}\color{red}{(n+1) a_{n+1}}x^{n}$$
The derivative:
$$\vec a'=D\vec a,\qquad D=\left( \begin{array}{ccccc} 0&0&0&0&\cdots\\ 1&0&0&0&\cdots\\ 0&2&0&0&\cdots\\ 0&0&3&0&\cdots\\ &&{\cdots}&& \end{array} \right)$$
Note: Monomials with high order are very bad for numerical calculation (better: use orthogonal basis)
Fourier basis $\displaystyle f(x)=\sum_{n=0}^Nc_ne^{inx}$
$$f'(x)=\sum_{n=0}^N\color{red} {inc_n}e^{inx},\quad\Rightarrow\quad D_{nm}=in\delta_{nm}$$
To switch between $c_n$ (Fourier) and $y_n$ (real) spaces uses Fast Fourier Transform (FFT). This matrix transformed to the real space is full not sparse.
Full matrices use all available information resulting in great accuracy $O(10^{-\alpha N})$.
def show_Fourier_coeffs(N = 30):
x = np.linspace(0, 2*np.pi, N, endpoint=False)
n = np.arange(N)
y = 1 + 2*np.sin(5*x)+0.5*np.cos(10*x) + np.exp(-11j*x)
y_fft = np.fft.fft(y)/N
plt.figure(figsize=(6,3), dpi=150)
plt.step(n+0.5, np.real(y_fft))
plt.step(n+0.5, np.imag(y_fft))
plt.grid(True)
plt.legend(['re$f$', 'im$f$'])
plt.xlabel('$n$')
plt.ylabel('$c_n$')
plt.show()
show_Fourier_coeffs(N = 30)
The vector of frequencies: $$\left[0,1, 2, \ldots, \frac{N}{2}, -\left(\frac{N}{2}-1\right), -\left(\frac{N}{2}-2\right),\ldots,-1\right]\frac{2\pi}{L} $$
Test for the function
$$f(x) = e^{\cos(x)}\;\;\text{for}\;\;x\in[0, 2\pi) $$$$f'(x) = -\sin(x)e^{\cos(x)}$$def f(x):
return np.exp(np.cos(x))
def df(x):
return -np.sin(x) * np.exp(np.cos(x))
def check_method(D, N):
x = np.linspace(0, 2*np.pi, N, endpoint=False)
h = x[1]-x[0]
F = f(x)
DF_acc = df(x)
return [N, D.__name__, D(F, DF_acc, h, N)]
def finite_diff(F, DF_acc, h, N):
n = np.arange(N)
DF = (F[n-2]-8*F[n-1]+8*F[(n+1)%N]-F[(n+2)%N]) / (12*h)
error = np.max(np.abs(DF-DF_acc))
return error
check_method(finite_diff, 100)
def spectral_diff(F, DF_acc, h, N):
g = np.empty([N,2], dtype=complex)
k = np.empty(N)
k[0:N//2+1] = np.arange(N//2+1)
n = np.arange(1,N//2)
k[-n] = -k[n]
ff = 1j*k*np.fft.fft(F)
DF = np.fft.ifft(ff)
error = np.max(np.abs(DF-DF_acc))
return error
#np.fft.ifft(1j*k*ff)/((N+1)/np.pi)
print(g)
N = 20
(check_method(finite_diff, N), check_method(spectral_diff, N)[1:])
def derivative_comparison():
K = 90
spec, fin, n = np.empty(K), np.empty(K), np.arange(K)
for k in n:
N = k+10
fin[k] = check_method(finite_diff, N)[2]
spec[k] = check_method(spectral_diff, N)[2]
plt.figure(figsize=(6,3), dpi=150)
plt.grid(True)
plt.semilogy(n+10, fin)
plt.semilogy(n+10, spec, 'r')
plt.legend(['finite difference', 'spectral method'])
plt.xlabel('$N$')
plt.ylabel('Error')
plt.show()
derivative_comparison()
from scipy.fftpack import dct, idct
def spectral_accuracy_problem():
N=1024;
x = np.arange(np.pi/(2*N), np.pi, np.pi/N) # (0,pi)
f = np.heaviside(x-0.25*np.pi, 0)-np.heaviside(x-0.75*np.pi, 0)
F=dct(f);
plt.figure(figsize=(10,3), dpi=120)
plt.subplot(1, 2, 1)
plt.grid(True)
n=np.arange(N)
plt.semilogy(n, F, 'ro');
plt.legend(['Fourier coefficients']);
F[50:]=0 # Truncate the Fourier series
f2=idct(F)/(2*N);
plt.subplot(1, 2, 2)
plt.grid(True)
plt.legend(['exact', 'truncated']);
plt.plot(x,f, label='exact', color="b");
plt.plot(x,f2, label='truncated', color='g');
plt.legend(['exact', 'truncated'])
plt.show()
spectral_accuracy_problem()
X = np.linspace(-1, 1, 11)
def l(n, x, Xn=X):
res = 1
for (k, xn) in enumerate(Xn):
if k!=n:
res *= (x-Xn[k])/(Xn[n]-Xn[k])
return res
def show_lagrange_example():
x = np.linspace(-1,1, 1000)
plt.grid(True)
plt.plot(x, l(1, x))
plt.plot(x, l(5, x))
plt.plot(x, l(7, x))
plt.legend(['$\ell_{1}(x)$', '$\ell_{5}(x)$', '$\ell_{7}(x)$'])
plt.plot(X, 0*X, 'o')
plt.plot([X[5], X[7], X[1]], [1, 1, 1], '+')
show_lagrange_example()
def show_interpolation_example(F, N=11):
x = np.linspace(-1,1, 1000)
X = np.linspace(-1,1, N)
Y = F(X)
y = F(x)
res = 0*x
for n in range(N):
res += l(n, x, X)*Y[n]
plt.grid(True)
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.ylim(-0.5, 1.5)
plt.plot(x, y)
plt.plot(x, res)
plt.plot(X, Y, "o")
plt.legend(['$f(x)$','$f_{aprx}(x)$','col. p.'])
plt.show()
return { 'max_error': np.max(np.abs(y-res))}
def f(x): return 1/(4*x**2+1)
plt.title("Lagrange approximation")
show_interpolation_example(f, N=11)
def f(x): return 1/(25*x**2+1)
plt.title("Lagrange interpolation")
show_interpolation_example(f, N=11)
def f(x): return 1/(25*x**2+1)
plt.title("Runge phenomenon")
show_interpolation_example(f, N=51)
For Lagrange interpolation large number of equidistant points lead to Runge phenomenon at the boundaries.
For fast vanishing functions it is possible to use sinc function approximation (no boundary)
$$\text{sinc}(x)=\frac{\sin \pi x}{\pi x}\;,\qquad \text{sinc}(n-m)=\delta_{nm}$$$$f(x)\approx \sum_{n=-\infty}^{\infty} f(x_n)\, \text{sinc}\left(\frac{x-x_n}{h}\right)\qquad x_n=nh$$def sinc(x): return np.sin(x*np.pi)/(x*np.pi)
import warnings
warnings.filterwarnings('ignore')
def show_sinc_interpolation(F, N):
x = np.linspace(-1,1, 1000)
X = np.linspace(-1,1, N)
h = X[1]-X[0]
Y = F(X)
y = F(x)
res = 0*x
for n in range(1,N-1):
res += Y[n] * sinc((x-X[n])/h)
plt.grid(True)
plt.xlabel("$x$")
plt.ylabel("$y$")
plt.ylim(-0.5, 1.5)
plt.plot(x, y)
plt.plot(x, res)
plt.plot(X, Y, "o")
plt.legend(['$f(x)$','$f_{aprx}(x)$','col. p.'])
plt.show()
return { 'max error': np.max(np.abs(y-res)), 'error without ends': np.max(np.abs(y-res)[100:-100])}
plt.title("Sinc interpolation")
show_sinc_interpolation(f, 11)
plt.title("Sinc interpolation")
show_sinc_interpolation(f, 51)
From the sinc interpolation a derivative can be calculated:
$$f'(x_m)\approx \sum_{n=-\infty}^{\infty} f_n\, \left.\frac{d}{dx}\text{sinc}\left(\frac{x-x_n}{h}\right)\right|_{x=x_m}=\sum_{n=-\infty}^{\infty} D_{mn}f_n$$$$D_{m-n}=\frac{1}{h}\left[\cdots, \frac13,-\frac12, 1, 0, -1, \frac12, -\frac13,\cdots\right]$$or
$$D_{mn} = \frac{(-1)^{m-n}}{m-n},\;m\neq n\;\;\text{and}\;\;D_{mm}=0$$Note: The matrix is full not sparse.
Other alternatives (More on that later??? - very important stuff):
cubic splines: piecewise defined approximation with cubic polynomial in each segment $Q_n(x), x\in[x_n,x_{n+1}]$ smoothly matched at grid points
$Q_{n-1}(x_n)=Q_n(x_n)$,
$Q'_{n-1}(x_n)=Q'_n(x_n)$,
$Q''_{n-1}(x_n)=Q''_n(x_n)$,
this requires solving linear algebraic equation with a tridiagonal matrix.