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!