# Why does SciPy eigsh() produce erroneous eigenvalues in case of harmonic oscillator?

I'm developing some larger code to perform eigenvalue computations of huge sparse matrices, in the context of computational physics. I test my routines against the simple harmonic oscillator in one dimension, since the eigenvalues are well known analytically. Doing so and comparing my own routines to SciPy's inbuilt solvers, I have come across the oddity displayed in the plot below. Here you can see the first 100 numerically computed eigenvalues $\lambda_{num}$ and analytical eigenvalues $\lambda_{ana}$

Around eigenvalue number 40, the numercial results start to diverge from the analytical ones. This doesn't surprise me (I won't go into why here, unless it comes up in the discussion). However, what is surprising to me is that eigsh() produces degenerate eigenvalues (around eigenvalue number 80). Why does eigsh() behave like that for even such a small number of eigenvalues?

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

#discretize x-axis
N = 100
xmin = -10.
xmax = 10.
accuracy = 1e-5
#stepsize
h = (xmax - xmin) / (N + 1.)
#exclude first and last points since we force wave function to be zero there
x = np.linspace(-10. + h,10. - h,N)
#create potential
V = x**2

def fivePoint(N,h,V):
C0 = (np.ones(N))*30. / (12. * h * h) + V
C1 = (np.ones(N)) * (-16.) / (12. * h * h)
C2 = (np.ones(N)) / (12. * h * h)
H = sp.spdiags([C2, C1, C0, C1, C2],[-2, -1, 0, 1, 2],N,N)
return H

H = myFunc.fivePoint(N,h,V)
eigval,eigvec = eigsh(H, k=N-1, which='SM', tol=accuracy)

#comparison analytical and numerical eigenvalues
xAxes = np.linspace(0,len(eigval)-1,len(eigval))
analyticalEigval = 2. * (xAxes + 0.5)
plt.figure()
plt.plot(xAxes,eigval, '+', label=r"$\lambda_{num}$")
plt.plot(xAxes,analyticalEigval, label=r"$\lambda_{ana}$")
plt.xlabel("Number of Eigenvalue")
plt.ylabel("Eigenvalue")
plt.legend(loc=4)
plt.title("eigsh()-method: Comparison of $\lambda_{num}$ and $\lambda_{ana}$")
plt.show()

• That's a very curious behavior. I will test it later today. Apr 14 '13 at 9:30
• I found the answer. In short: my thinking was wrong. The analytical solutions of the harmonic oscillator (HOSZ) are valid without any spacial restrictions. However, in the code above, my box runs from -10 to 10, so this puts a boundary condition on the numerical solutions. Consequently, eigsh() solves the system it is given correctly. At around n = 50 (with n being the principal Quantum number), the analytical solutions do not fit inside the -10, 10 box anymore. Now (after some thinking), this seems obvious. However, I didn't see that while building and testing the code.
– seb
Apr 14 '13 at 10:51
• this still doesn't explain the degeneracy though, does it?
– seb
Apr 14 '13 at 10:54

Nonetheless, no method is perfect, especially if you're trying to compute the entire spectrum of your matrix. So if you're trying to get the 50 smallest eigenvalues, you may be better off approximating the wave function by a vector with 100 elements and only asking eigsh() for the first 50 energy levels, rather than using a vector with 50 points and asking for all of the eigenvalues.