14
$\begingroup$

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?

enter image description here

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()
$\endgroup$
  • $\begingroup$ That's a very curious behavior. I will test it later today. $\endgroup$ – Rafael Reiter Apr 14 '13 at 9:30
  • $\begingroup$ 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. $\endgroup$ – seb Apr 14 '13 at 10:51
  • $\begingroup$ this still doesn't explain the degeneracy though, does it? $\endgroup$ – seb Apr 14 '13 at 10:54
11
$\begingroup$

The degeneracy of some eigenvalues looks to me like the hallmark of the breakdown of the Lanczos algorithm. The Lanczos algorithm is one of the more commonly used methods to approximate the eigenvalues and eigenvectors of Hermitian matrices; it's what scipy.eigsh() uses, through a call to the ARPACK library.

In exact arithmetic, the Lanczos algorithm produces a set of orthogonal vectors, but in floating point arithmetic these can fail to be orthogonal and even become linearly dependent. The really annoying thing is that this loss of orthogonality happens precisely when one of the approximate eigenvalues has converged to one of the real eigenvalues -- the algorithm sabotages itself, so to speak. The result is that you'll get some spurious pairs of nearby eigenvalues. There are various fixes for this, for example using Gram-Schmidt to force any converged eigenvectors to be orthogonal at every step.

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.

If you want to read more, look at Yousef Saad's Numerical Methods for Large Eigenvalue Problems.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.