氢原子的Schrodinger方程：为什么numpy在显示scipy时显示错误的解？

Question

我编写了一段代码来解决一维Schrodinger方程。尽管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是否包含一个例程，我可以安全地使用该例程来解决特征值问题，而不论其潜力如何？

提前致谢！

Answer 1

which='SM'in中的参数eigsh告诉函数找到幅度k最小的特征值。最大负特征值约为-15.22，但不会包含在的结果中，因为还有许多其他正负特征值，其大小小于15.22。如果您想将结果与（或更好）进行比较，请使用。这告诉我们找到代数最小值。如果这样做，则您使用进行计算的将会与所计算出的前6个特征值一致。eigshnumpy.linalg.eignumpy.linalg.eighwhich='SA'eigshvalseigshnumpy.linalg.eigh

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

GS = vecs[:,1]

使图的精确度和数值结果一致。