Finite Differences

In [1]:
import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as pt

Part 1: Examining the Differentiation Matrix

In [2]:
degree = 2
h = 0.25

# Assume even degree so that there's a well-defined middle node.
assert degree % 2 == 0

nodes = np.linspace(-h/2, h/2, degree+1)
nodes
Out[2]:
array([-0.125,  0.   ,  0.125])

Now construct V (the generalized Vandermonde) and Vprime (the generalized Vandermonde for the derivatives):

In [3]:
V = np.array([
    nodes**i
    for i in range(degree+1)
]).T
In [4]:
def monomial_deriv(i, x):
    if i == 0:
        return 0*x
    else:
        return i*nodes**(i-1)

Vprime = np.array([
    monomial_deriv(i, nodes)
    for i in range(degree+1)
]).T

Combine them to form the derivative matrix:

In [6]:
diff_mat = Vprime.dot(la.inv(V))
diff_mat
Out[6]:
array([[-12.,  16.,  -4.],
       [ -4.,   0.,   4.],
       [  4., -16.,  12.]])

Let's say we only care about the derivative at the middle node:

In [33]:
finite_difference_weights = diff_mat[degree//2]
finite_difference_weights
Out[33]:
array([-4.,  0.,  4.])
  • What have we learned?
  • What formula does this amount to?
  • How do these weights change if we change $h$?
  • What formula does this amount to, really?
  • What happens if we change the degree?
  • What happens if we shift all nodes?
In [10]:
# * We could have left the middle point out. :)
# * -4*f(x-0.25) + 4*f(x+0.25)
# * They scale with 1/h, as you might expect.
# * (f(x-h/2) + f(x+h/2))/h
# * We get a more complicated (but more accurate) formula (with more source nodes)
# * The weights remain the same.

Part 2: Using finite difference formulas

In [55]:
def f(x):
    return np.sin(4*x)
def df(x):
    return 4*np.cos(4*x)
In [56]:
x = np.arange(10) * 0.125
pt.plot(x, f(x), "o-")
Out[56]:
[<matplotlib.lines.Line2D at 0x7f761043a198>]

Now use the weights to compute the finite difference derivative as deriv:

In [57]:
fdw = finite_difference_weights

fx = f(x)
deriv = np.zeros(len(x)-2)
for i in range(1, 1+len(deriv)):
    deriv[i-1] = fx[i-1]*fdw[0] + fx[i]*fdw[1] + fx[i+1]*fdw[2]

Now plot the finite difference derivative:

In [59]:
pt.plot(x[1:-1], df(x[1:-1]), label="true")
pt.plot(x[1:-1], deriv, label="FD")
pt.legend()
Out[59]:
<matplotlib.legend.Legend at 0x7f7610378860>