# Correct eigenfunctions of Laplace operator by Finite Differences

I am trying to compute the eigenfunctions of the Laplace operator, i.e. finding $u$ in

$$-\nabla^2 u = \lambda u .$$

For now I am trying to do this in 1D, so

$$\nabla^2 = \partial_{xx} .$$

I am using the finite difference method. My boundary conditions are $u=0$ at $x = -1, 1$. Here is my code, and the results:

import numpy as np
from scipy.sparse.linalg import eigs
import matplotlib.pyplot as plt

n = 200
h = 2/(n-1)     # domain for x and y is [-1, 1]

L = np.diag(np.ones(n-1), k=-1) + np.diag(-2*np.ones(n)) + \
np.diag(np.ones(n-1), k=1)
L *= -1/h**2

eigvals, eigvecs = eigs(L)

eig = np.real(eigvecs[:, 0])

x = np.linspace(-1, 1, num=n)

plt.plot(x, eig)
plt.show()


Can anyone see what I am doing wrong, or give suggestions as to how I might be able to fix my code?

You should specify the eigenvalues you want with which="SM", for example.

Check the following snippet. I also changed the solver, since your system is symmetric.

import numpy as np
from scipy.sparse.linalg import eigsh
import matplotlib.pyplot as plt

n = 200
h = 2/(n-1)     # domain for x and y is [-1, 1]

L = np.diag(np.ones(n-1), k=-1) - np.diag(2*np.ones(n)) + \
np.diag(np.ones(n-1), k=1)
L = -L/h**2

eigvals, eigvecs = eigsh(L, which="SM")

eig = eigvecs[:, 0]

x = np.linspace(-1, 1, num=n)

plt.plot(x, eig)
plt.show()


• Thanks for your answer - although I do have a follow-up question. If I don't specify the smallest eigenvalues, then the number of "spikes" in the solution increases with n, with every point alternating positive and negative, rather than seeming to converge towards a single eigenvector. Is this consistent with the explanation that the solver is finding higher order solutions? – islanss Jun 16 '17 at 20:06
• I think that is the expected result. Your discrete version of the operator approximates better smaller eigenvalues/eigenvectors. I think that @SpencerBryngelson answer makes sense there. – nicoguaro Jun 16 '17 at 21:08