氢原子的薛定谔方程：为什么numpy显示出错误的解，而scipy没有？

Question

我写了一段代码来解一维薛定谔方程。而纽比.利纳格.艾格（）对于谐振子，常规方法已经很好地工作了，它似乎为库仑势增加了一个假解。另一方面，西皮的稀疏.linalg.eigsh（）似乎做得很好。这是我的剧本：

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

N = 500
x0 = 8
xMin, xMax = -x0, x0
xstep = (xMax - xMin) / (N - 1)
x = np.linspace(xMin, xMax, N)

k = np.array([np.ones(N-1),-2*np.ones(N),np.ones(N-1)])
offset = [-1,0,1]
Lapl = diags(k,offset).toarray() / (xstep**2)
T = -Lapl *.5

V = -1./(np.abs(x)) #Coulomb
#V = .5 * x**2 #HO

H = T.copy()
for i in range(N):
    H[i,i] += V[i]

#vals, vecs = np.linalg.eig(H)
vals, vecs = eigsh(H, which='SM')
idx = vals.argsort()[::1]
vals, vecs = vals[idx], vecs[:,idx]

#exact solution
Hatom = (np.pi)**(-1./2) * np.exp(- np.abs(x)) * np.abs(x) * np.sqrt(4 * np.pi)
norm = np.sum(Hatom**2)
Hatom = Hatom / np.sqrt(norm)

#numerical solution
GS = vecs[:,0] #0th is the gs if using sp's eigsh, but is spurious if using np.linalg.eig
normGS = np.sum(GS**2)
GS = GS / np.sqrt(normGS)

plt.plot(x, Hatom**2, label = 'exact')
plt.plot(x, GS**2, label = 'numeric')
plt.legend()
plt.show()

print( np.round(vals[:10], 4) )

这就产生了下面的图（我也有困难直接在这里显示图片，很抱歉！）公司名称：

• 使用numpy:np_0th_eigenvector_spurious &；np_1th_eigenvector_gs

• 使用scipy:sp_gs

我希望这来自于对库仑势奇异性的不同处理（尽管我选择了偶数个点来避免x=0），因为numpy和scipy对于谐振子和Morse势都很好（为了简洁起见，这里不重复）。然而，这使得如何处理任意电位变得很棘手！在

此外，库仑势场的本征值与1/n^2序列相距甚远（最小值来自于使用numpy）：

^{pr2}$

我做错什么了？numpy/scipy是否包含一个我可以安全地用来解决特征值问题的例程，而不管其潜力如何？在

提前谢谢！在

Answer 1

在^{}中的自变量which='SM'告诉函数以最小的大小找到k特征值。最负的特征值约为-15.22，但这不会包括在eigsh的结果中，因为还有许多其他正负本征值小于15.22。如果要将结果与^{}（或更好的，^{}）的结果进行比较，请使用which='SA'。它告诉eigsh找到代数上最小的值。如果您这样做，那么您用eigsh计算的vals将与numpy.linalg.eigh计算的前6个特征值一致。在

在进行此更改后，您必须更改要绘制到的数值计算特征向量的选择

GS = vecs[:,1]

使精确结果和数值结果的曲线一致。在

Answer 2

按照@Warren所说的whichk特征向量和特征值，你也可以利用矩阵H的三对角结构，这里使用scipy.linalg.eigh_tridiagonal，而不是{}，这给你：

from scipy.linalg import eigh_tridiagonal
d = np.diag(H)
e = np.diag(H,k=1)
vals, vecs = eigh_tridiagonal(d,e)
idx = vals.argsort()[::1]
vals, vecs = vals[idx], vecs[:,idx]

并相应地选择GS = vecs[:,1]以获得相同的结果，计算速度提高2倍：

^{pr2}$