Ch19.5 Part 3: Numerical Differentiation

I just went to an emotional wedding.

Even the cake was in tiers.

Standard numerical methods for computing derivatives correspond to differentiating interpolating polynomials for values of \( f(x) \) on nodes \( x_k \).

• Two-point rules uses linear interpolation on nodes.
• Three-point rules uses quadratic interpolation on nodes.

For two nodes \( \{ x_0, x_1 \} \), where \( x_1 = x_0 + h \), the Lagrange interpolating polynomial and its derivative are

\[ \small{ \begin{aligned} P(x) & = f(x_0)\frac{x - x_1}{x_0-x_1} + f(x_1) \frac{x - x_0}{x_1 - x_0} \\ & = f(x_0)\frac{x - x_0 - h}{- h} + f(x_0 + h) \frac{x - x_0}{h} \\ P'(x_0) & = f(x_0)\frac{1}{- h} + f(x_0 + h) \frac{1}{h} \\ & = \frac{f(x_0 + h)-f(x_0)}{h} \end{aligned} } \]

• Two-Point Formula
• This formula looks both ahead and behind:

\[ f'(x) \cong \frac{f(x+h)-f(x-h)}{2h} \]

Forward:

\[ f'(x) \cong \frac{-3f(x)+4f(x+h)-f(x+2h)}{2h} \]

Backward:

\[ f'(x) \cong \frac{3f(x)-4f(x-h)+f(x-2h)}{2h} \]

Centered:

\[ f'(x) \cong \frac{f(x+h)+f(x-h)}{2h} \]

For three nodes \( x_0, x_1 = x_0 + h, x_2 = x_0 + 2h \), the Lagrange interpolating polynomial and its derivative are

\[ \scriptsize{ \begin{aligned} P(x) & = f(x_0)\frac{(x - x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} + f(x_1) \frac{(x - x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} + f(x_2)\frac{(x - x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} \\ & = f(x_0)\frac{(x - x_0-h)(x-x_0-2h)}{(-h)(-2h)} + f(x_0 +h) \frac{(x - x_0)(x-x_0-2h)}{(h)(-h)} \\ & \hspace{1in} + f(x_0 + 2h)\frac{(x - x_0)(x-x_0-h)}{(2h)(h)} \\ \\ P'(x_0) & = \mathrm{Product \, rule, \, algebra, \, plug \, in \, } x = x_0 \\ \\ & = \frac{-3f(x_0)+4f(x_0+h)-f(x_0+2h)}{2h} \end{aligned} } \]

def fwd(f,a,h=0.01):
    fp = (f(a+h)-f(a))/h
    return fp

import numpy as np
f = np.sin
fwd(f,0)

0.9999833334166665

fwd(f,0,-0.0001) #bwd diff

0.9999999983333334

def fwd(f,a,h=0.01):
    fp = (f(a+h)-f(a))/h
    return fp

import numpy as np
f = np.sin
fwd(f,np.pi/2,0.0001)

-4.999999969612645e-05

fwd(f,np.pi/2,-0.1) #bwd diff

0.049958347219742905

def bwd3(f,a,h=0.0001):
    fp = (3*f(a)-4*f(a-h)+f(a-2*h))/(2*h)
    return fp

import numpy as np
f = np.sin
x = np.array([0,np.pi/2])
bwd3(f,x)

array([1., 0.])

def sym(f,a,h=0.0001):
    fp = (f(a+h)-f(a-h))/(2*h)
    return fp

import numpy as np
f = np.sin
x = np.array([0,np.pi/2])
sym(f,x)

array([1., 0.])

• Suppose \( f(x) \) is given in terms of input and output data vectors rather than as an analytic expression.
• Then \( f'(x) \) will consist of input and output vectors whose lengths are one less than originals (forward difference) or two less (centered) because of boundary constraints.

\[ \small{ \begin{aligned} f'(x) &\cong \frac{f(x+h)-f(x)}{h} \\ f'(x) &\cong \frac{f(x+h)-f(x-h)}{2h} \end{aligned} } \]

xdata = [1, 1.1, 1.2, 1.3, 1.4, 1.5]
ydata = [2, 3, 4, 4.5, 5, 6]

• It can be easier sometimes to use np.arange(a,b,h) instead of np.linspace(a,b,n), because you specify step size \( h \) rather than number of subintervals \( n \).
  
• Also, see command below for omitting last entry in a vector.

import numpy as np
x = np.arange(0,1,0.1)
print(x)
print(x[:-1])

[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]

[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8]

For this example, we use the np.diff command to evaluate the forward differences of \( \sin(2 \pi x) \) on a vector of nodes.

import numpy as np
import matplotlib.pyplot as plt
h = 0.01
x = np.arange(0,2*np.pi,h)
f = np.sin(x)
fp = np.diff(f)/h
plt.figure(figsize=(10,5))
plt.plot(x,f,'b',label='y = f(x)')
plt.plot(x[:-1],fp,'r',label="y = f'(x)")
plt.title("f(x) = sin(x) in blue, and f'(x) = cos(x) in red")
plt.legend(loc='best')
plt.show()

import numpy as np
from scipy.interpolate import interp1d

#Command for cubic spline spline polynomial p(x)
f = interp1d(xdata, ydata, kind = 'cubic')
N = 100
x = np.linspace(xdata[0],xdata[n-1],N)
y = f(x)

h = 0.01
fp = np.zeros(N)
fp[0] = fwd(f,x[0])        #fwd diff
fp[N-1] = fwd(f,x[N-1],-h) #bwd diff
for k in range(1,N-1):
    fp[k] = sym(f,x[k],h)  #ctr diff

import numpy as np
from scipy.interpolate import interp1d

#Command for cubic spline spline polynomial p(x)
f = interp1d(xdata, ydata, kind = 'cubic')
N = 100
x = np.linspace(xdata[0],xdata[n-1],N)
y = f(x)

h = 0.01
fp = np.zeros(N)
fp[0] = fwd(f,x[0])            #fwd diff
fp[N-1] = fwd(f,x[N-1],-h)     #bwd diff
fp[1:N-1] = sym(f,x[1:N-1],h)  #ctr diff

Find error bound for forward (backward) difference.

\[ \small{ f(x+h) = f(x) + f'(x)\left((x+h)-x)\right) + \frac{f''(x)}{2}\left((x+h)-x)\right)^2 + \ldots } \]

It follows that

\[ \small{ \begin{aligned} f(x+h) - f(x) & = f'(x)h + \frac{f''(x)}{2}h^2 + \ldots \\ f'(x) & = \frac{f(x+h) - f(x)}{h} - \frac{f''(x)}{2}h + \ldots \\ f'(x) & = \frac{f(x+h) - f(x)}{h} - \frac{f''(c(x))}{2}h, \,\, \,\, c(x) \in [x,x+h] \\ \end{aligned} } \]

Taylor series again, forward and backward: \[ \small{ \begin{aligned} f(x+h) & = f(x) + f'(x)\left((x+h)-x)\right) + \frac{f''(x)}{2}\left((x+h)-x)\right)^2 + \ldots \\ &= f(x) + f'(x)h + \frac{f''(x)}{2}h^2 + \frac{f'''(x)}{6}h^3 + \ldots \\ f(x-h) & = f(x) - f'(x)h + \frac{f''(x)}{2}h^2 - \frac{f'''(x)}{6}h^3 + \ldots \\ \end{aligned} } \]

It follows that

\[ \small{ \begin{aligned} f'(x) & = \frac{f(x+h) - f(x-h)}{2h} - \frac{f'''(x)}{6}h^3 + \ldots \\ f'(x) & = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(c(x)), \,\, \,\, c(x) \in [x-h,x+h] \\ \end{aligned} } \]

Here are the forward and centered difference formulas and their error terms in the 2-point case:

\[ \small{ \begin{aligned} f'(x) & = \frac{f(x+h) - f(x)}{h} - \frac{f''(c(x))}{2}h, \,\, \,\, c(x) \in [x,x+h] \\ \\ f'(x) & = \frac{f(x+h) - f(x-h)}{2h} - \frac{h^2}{6}f'''(c(x)), \,\, \,\, c(x) \in [x-h,x+h] \\ \end{aligned} } \]

• We can expect the centered difference formula to be more accurate.
• The forward and backward difference formulas can be used at the endpoints.