我写了一段代码来解决一维薛定谔方程。虽然numpy.linalg.eig((例程对于谐波振荡器工作正常，但它似乎为库仑势增加了一个杂散解。另一方面，Scipy的sparse.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 对于谐波振荡器和莫尔斯势都很好(为了简洁起见，这里不复制(。然而，这使得如何处理任意电位变得棘手！
更重要的是，库仑势地的特征值与 1/n^2 序列相去甚远(最低的特征值来自使用 numpy(：
vals: [-15.220171 -0.500363 -0.331042 -0.085621 -0.02475 0.242598
0.344308 0.741555 0.885751 1.402606]
我在这里做错了什么？numpy/scipy 是否包含我可以安全地用于解决特征值问题的例程，无论潜力如何？
提前感谢！
eigsh中which='SM'的参数告诉函数找到具有最小量级的k特征值。 最负的特征值大约是-15.22，但这不会包含在eigsh的结果中，因为还有许多其他正特征值和负特征值的大小小于15.22。 如果要将结果与numpy.linalg.eig(或更好的numpy.linalg.eigh(的结果进行比较，请使用which='SA'。 这告诉eigsh找到代数最小的值。 如果这样做，则使用eigsh计算的vals将与numpy.linalg.eigh计算的前 6 个特征值一致。
进行此更改后，您必须更改要绘制的数值计算特征向量的选择。
GS = vecs[:,1]
使精确结果和数值结果的图一致。
按照@Warren所说的whichk 个特征向量和特征值来查找，您还可以利用矩阵H的三对角线结构，通过在此处使用 scipy.linalg.eigh_tridiagonal 而不仅仅是np.linalg.eigh为您提供：
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 倍：
%timeit scipy.linalg.eigh_tridiagonal(d,e)
>>> 10.7 ms ± 37.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit numpy.linalg.eigh(H)
>>> 19.4 ms ± 28.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit scipy.sparse.linalg.eigsh(H, which='SA')
>>> 115 ms ± 1.98 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)