I am trying to solve for the eigenfunctions of a (1D) differential operator using finite differences:
$$A \, f(x) = \lambda f(x)$$
Here is an example in Python where $A = \partial_x^4$:
import numpy as np import matplotlib.pyplot as plt from scipy.sparse import spdiags from scipy import linalg def derivative_matrix(n): stencil = np.array((1,-4,6,-4,1)) diags = range(-2,3) bands = np.tile(stencil,(n,1)).T A = spdiags(bands,diags,n,n).todense() A = np.array(A,dtype=float) return A if __name__=="__main__": n = 100 A = derivative_matrix(n) vals,vecs = linalg.eig(A) plt.plot(np.real(vecs[:,-3:None])) plt.show()
This gives me really nice eigenfunctions:
My questions is: what do I have to change in the implementation to change the boundary condition of $f$? Currently, the method seems to implicitely assume I want $f(x)=f'(x)=0$. But what if I instead wanted $f'(x)=f''(x)=0$? How would I implement that?
Any help will be much appreciated. Thanks in advance!