# coding: utf-8 # # 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 # 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 # 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 # * 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-") # 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()