Dlaczego SciPy eigsh () wytwarza błędne wartości własne w przypadku oscylatora harmonicznego?

Rozwijam jakiś większy kod do wykonywania obliczeń wartości własnych ogromnych rzadkich macierzy w kontekście fizyki obliczeniowej. Moje procedury sprawdzam na prostym oscylatorze harmonicznym w jednym wymiarze, ponieważ wartości własne są dobrze znane analitycznie. Robiąc to i porównując własne procedury z wbudowanymi rozwiązaniami SciPy, natknąłem się na osobliwość pokazaną na poniższym wykresie. Tutaj możesz zobaczyć pierwsze 100 obliczonych liczbowo wartości własnych i analitycznych wartości własnych λ a n $\lambda_{num}$$\lambda_{ana}$

Około wartości własnej 40, wyniki liczbowe zaczynają odbiegać od wyników analitycznych. Nie zaskakuje mnie to (nie będę się tutaj zastanawiać, chyba że pojawi się w dyskusji). Jednak zaskakuje mnie to, że eigsh () wytwarza zdegenerowane wartości własne (wokół wartości własnej 80). Dlaczego eigsh () zachowuje się tak nawet w przypadku tak małej liczby wartości własnych?

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()

Degeneracja niektórych wartości własnych wydaje mi się znakiem rozpoznawczym załamania się algorytmu Lanczosa . Algorytm Lanczosa jest jedną z najczęściej używanych metod przybliżania wartości własnych i wektorów własnych macierzy hermitowskich; używa tego scipy.eigsh () poprzez wywołanie biblioteki ARPACK .

Dokładnie w arytmetyce algorytm Lanczosa tworzy zbiór wektorów ortogonalnych, ale w arytmetyki zmiennoprzecinkowej mogą one nie być ortogonalne, a nawet stać się zależne liniowo. Naprawdę denerwujące jest to, że utrata ortogonalności następuje dokładnie wtedy, gdy jedna z przybliżonych wartości własnych zbiegnie się z jedną z rzeczywistych wartości własnych - że tak powiem, algorytm sabotuje sam. W rezultacie otrzymasz kilka fałszywych par pobliskich wartości własnych. Istnieją różne poprawki tego, na przykład użycie Grama-Schmidta do wymuszenia, aby każdy zbieżny wektor własny był prostopadły na każdym kroku.

Niemniej jednak żadna metoda nie jest idealna, szczególnie jeśli próbujesz obliczyć całe spektrum macierzy . Więc jeśli próbujesz uzyskać 50 najmniejszych wartości własnych, możesz lepiej zbliżyć funkcję fali do wektora o 100 elementach i poprosić tylko o eigsh()pierwsze 50 poziomów energii, zamiast używać wektora o 50 punktach i poprosić o wszystko wartości własnych.

Jeśli chcesz dowiedzieć się więcej, spójrz na numeryczne metody Yousefa Saada dla problemów z dużymi wartościami własnymi .