python - 氢原子的薛定谔方程 : why is numpy displaying a wrong solution while scipy isn't?

Question

我编写了一段代码来求解一维薛定谔方程。虽然 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

eigsh 中的参数 which='SM'告诉函数找到 k 个具有最小 magnitude 的特征值。最负的特征值约为 -15.22，但这不会包含在 eigsh 的结果中，因为还有许多其他正负特征值的幅度小于 15.22。如果您想将结果与 numpy.linalg.eig 的结果进行比较(或更好，numpy.linalg.eigh )，使用 which='SA'。这告诉 eigsh 找到代数上最小的值。如果你这样做，那么你用 eigsh 计算的 vals 将与 numpy.linalg.eigh 计算的前 6 个特征值一致。

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

GS = vecs[:,1]

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