Почему SciPy eigsh () выдает ошибочные собственные значения в случае гармонического осциллятора?

15

Я разрабатываю некоторый более крупный код для выполнения вычислений по собственным значениям огромных разреженных матриц в контексте вычислительной физики. Я проверяю свои процедуры против простого гармонического осциллятора в одном измерении, поскольку собственные значения хорошо известны аналитически. Делая это и сравнивая мои собственные процедуры со встроенными решателями SciPy, я столкнулся со странностью, показанной на графике ниже. Здесь вы можете увидеть первые 100 численно вычисленных собственных значений и аналитических собственных значений λ a n aλNUмλaNa

Вокруг собственного значения 40 числовые результаты начинают расходиться с аналитическими. Меня это не удивляет (я не буду вдаваться в причину, если только это не будет обсуждаться). Однако, что удивляет меня, так это то, что eigsh () создает вырожденные собственные значения (около 80 собственных значений). Почему eigsh () ведет себя так даже для такого небольшого числа собственных значений?

введите описание изображения здесь

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()
СЕБ
источник
Это очень любопытное поведение. Я проверю это позже сегодня.
Рафаэль Райтер
Я нашел ответ. Короче говоря, мое мышление было неверным. Аналитические решения гармонического осциллятора (HOSZ) действительны без каких-либо пространственных ограничений. Тем не менее, в приведенном выше коде мой ящик работает от -10 до 10, так что это ставит граничное условие для численных решений. Следовательно, eigsh () решает, какая система задана правильно. При n = 50 (при этом n является основным квантовым числом), аналитические решения больше не помещаются в поле -10, 10. Теперь (после некоторых размышлений) это кажется очевидным. Однако я не видел этого при сборке и тестировании кода.
СЕБ
это все еще не объясняет вырождение, не так ли?
СЕБ

Ответы:

12

Вырождение некоторых собственных значений выглядит для меня как отличительная черта нарушения алгоритма Ланцоша . Алгоритм Ланцоша является одним из наиболее часто используемых методов аппроксимации собственных значений и собственных векторов эрмитовых матриц; это то, что scipy.eigsh () использует, через вызов библиотеки ARPACK .

В точной арифметике алгоритм Ланцоша производит набор ортогональных векторов, но в арифметике с плавающей запятой они могут не быть ортогональными и даже стать линейно зависимыми. Действительно раздражает то, что эта потеря ортогональности происходит именно тогда, когда одно из приближенных собственных значений приблизилось к одному из реальных собственных значений - алгоритм, так сказать, саботирует сам себя. В результате вы получите несколько ложных пар близлежащих собственных значений. Для этого есть различные способы, например, использование Грам-Шмидта, чтобы заставить любые сходящиеся собственные векторы быть ортогональными на каждом шаге.

Тем не менее, ни один метод не является идеальным, особенно если вы пытаетесь вычислить весь спектр вашей матрицы . Поэтому, если вы пытаетесь получить 50 наименьших собственных значений, вам может быть лучше аппроксимировать волновую функцию вектором с 100 элементами и запрашивать только eigsh()первые 50 уровней энергии, а не использовать вектор с 50 точками и запрашивать все из собственных значений.

Если вы хотите узнать больше, посмотрите на Численные методы Юсефа Саада для больших задач на собственные значения .

Даниэль Шаперо
источник