现在在自学赖文的量子化学，最近看到简谐振子这一部分。书中以简谐振子为例，介绍了使用Numerov方法数值求解Schroedinger方程的方法，并提供了一个C++的代码来实现。因为我不会C++，所以就使用正在学的Python来尝试写一个（刚学到函数）。

原理

我在维基上搜了一下，Numerov方法不只用于量子力学。但是涉及的其他的领域我不了解，所以无法介绍。总之就是单纯的讲讲求解一维谐振子波函数的方法。

$$\hat{H}=\hat{T}+\hat{V}=-\frac{{\hslash }^2}{2m}\frac{d^2}{d{x}^2}+\frac{1}{2}k{x}^2=-\frac{{\hslash }^2}{2m}\frac{d^2}{d{x}^2}+2{\pi }^2{\nu }^{2}m{x}^2=-\frac{{\hslash }^2}{2m}(\frac{d^2}{d{x}^2}-{\alpha }^2{x}^2)$$ $$\alpha =\frac{2\pi \nu m}{\hslash }$$

一维量子谐振子的Schroedinger方程： $$\frac{d^2\psi }{d{x}^2}+(2mE{\hslash }^{-2}-{\alpha }^2{x}^2)\psi =0$$

一维谐振子的解析解的求解不算很难，这里也不做介绍了。

Numerov方法是基于递推，将我们所探讨的区域划分为多个小的间隔$s$。若我们取最左端的点的坐标为$x_0$，向右的点分别为$x_1=x_0+s,x_2=x_1+s,\dots ,x_{\text{max}}$。设${\psi }_n,{\psi }_{n+1},{\psi }_{n-1}$分别表示$\psi$在$x_n,x_n+s,x_n-s$处的值：

$$\psi (x_n)={\psi }_n,\psi (x_n+s)={\psi }_{n+1},\psi (x_n-s)={\psi }_{n-1}$$

将Schroedinger方程写为：

$${\psi }^{″}=G\psi ,G=m{\hslash }^{-2}[2V(x)-E]$$

把$\psi (x_n+s)$和$\psi (x_n-s)$进行Taylor展开：

$$\psi (x_n+s)=\psi (n)+{\psi }^{\prime }(n)\cdot s+\frac{{\psi }^{″}(n)}{2}\cdot s^2+\frac{{\psi }^{‴}(n)}{3!}\cdot s^3+\frac{{\psi }^{\text{(4)}}(n)}{4!}\cdot s^4+\frac{{\psi }^{\text{(5)}}(n)}{5!}\cdot s^5+\dots$$

$$\psi (x_n-s)=\psi (n)-{\psi }^{\prime }(n)\cdot s+\frac{{\psi }^{″}(n)}{2}\cdot s^2-\frac{{\psi }^{‴}(n)}{3!}\cdot s^3+\frac{{\psi }^{\text{(4)}}(n)}{4!}\cdot s^4-\frac{{\psi }^{\text{(5)}}(n)}{5!}\cdot s^5+\dots$$

将两式相加（忽略后面的项）：

$${\psi }_{n+1}=-{\psi }_{n-1}+2{\psi }_n+{\psi }_n^{″}\cdot s^2+\frac{{\psi }_n^{\text{(4)}}}{12}\cdot s^2$$

用${\psi }^{″}$代替$\psi$，并两边乘以$s^2$：

$${\psi }_{n+1}^{″}\cdot s^2=-{\psi }_{n-1}^{″}\cdot s^2+2{\psi }_n^{″}\cdot s^2+{\psi }_n^{\text{(4)}}\cdot s^4+\frac{{\psi }_n^{\text{(6)}}}{12}\cdot s^6$$

忽略含6阶导的项：

$${\psi }_n^{\text{(4)}}\cdot s^4={\psi }_{n+1}^{″}\cdot s^2+{\psi }_{n-1}^{″}\cdot s^2-2{\psi }_n^{″}\cdot s^2$$

将上式代入之前的式子，并代入${\psi }^{″}=G\psi$：

$${\psi }_{n+1}=-{\psi }_{n-1}+2{\psi }_n+G_n{\psi }_n{s}^2+\frac{1}{12}[G_{n+1}{\psi }_{n+1}{s}^2+G_{n-1}{\psi }_{n-1}{s}^2-2G_n{\psi }_n{s}^2]$$

所以：

$${\psi }_{n+1}=\frac{2{\psi }_n-{\psi }_{n-1}+5G_n{\psi }_n{s}^2/6+G_{n-1}{\psi }_{n-1}{s}^2/12}{1-G_{n+1}{s}^2/12}$$ $$G_n=G(x_n)=m{\hslash }^{-2}[2V(x_n)-2E]$$

由此我们得到了该递推关系式。为了用该式求得波函数，我们首先需要得到${\psi }_0,{\psi }_1$。 根据经典力学，$E<V$的区域为禁阻区域；而在量子力学中，由于隧穿效应，在$E<V$的区域，波函数$\psi$并不为零。但是进入经典禁阻区域越深，$\psi$越趋于零。因此我们假定$\psi (x_0)={\psi }_0=0$。而对于$x_1$点，假设波函数在该点为一很小的值，比如可以取0.0001。事实上取0.001也是可以的，这相当于将波函数乘以10，通过归一化处理后并没有什么区别。

要进行递推，显然我们还需要猜测一个本征能量值$E_{\text{guess}}$。经过递推最终得到$\psi (x_{\text{max}})$。根据品优波函数平方可积或对称性的要求，$\psi (x_{\text{max}})$应当趋近于零。若$\psi (x_{\text{max}})$远远偏离0，说明我们猜测的能量本征值不够合理，需要重新进行猜测。

有了上述的方法，我们尝试用程序进行实现。

程序实现

为了方便地进行计算，首先我们要做的一件事是将前面涉及的式子去量纲化（更准确的说法是将其转换为自然单位制）。

具体的变换不写了，并不算难。最终变换得到的Schroedinger方程为： $$\frac{d^2{\psi }_r}{d{x}_r^2}=(2V_r-2E_r){\psi }_r$$

$${\psi }_r^{″}=G_r{\psi }_r,G_r=2V_r-2E_r$$

上式中：$V_r=\frac{1}{2}{x}_r^2$。

经过处理，这些式子中的各个量的量纲均为1。

然后就是代码实现了。

（瞎画了一个流程图）

然后是python代码：

Other
开启自动换行复制
X={"x0":0.0}
P={"p0":0.0}
G={"g0":0.0}
nn=0

x=float(input("Enter initial xr:"))
s=float(input("Enter the increment sr:"))
m=int(input("Enter the number of intervals m:"))
E=float(input("Enter the reduced energy Er:"))

for i in list(range(m+1)):
	X["x%s" % i]=0.0 #储存x

for i in list(range(m+1)):
	P["p%s" % i]=0.0 #储存psi
	
for i in list(range(m+1)):
	G["g%s" % i]=0.0 #储存g

P["p1"]=0.0001
X["x0"]=x
X["x1"]=X["x0"]+s
G["g0"]=X["x0"]*X["x0"]-2*E
G["g1"]=X["x1"]*X["x1"]-2*E
ss=s*s/12

for i in range(1,m+1):
	X["x%s" % (i+1)]=X["x%s" % i]+s
	G["g%s" % (i+1)]=X["x%s" % (i+1)]*X["x%s" % (i+1)]-2*E
	P["p%s" % (i+1)]=(-P["p%s" % (i-1)]+2*P["p%s" % i]+10*G["g%s" % i]*P["p%s" % i]*ss+G["g%s" % (i-1)]*P["p%s" % (i-1)]*ss)/(1-G["g%s" % (i+1)]*ss)
	if (P["p%s" % i]*P["p%s" % (i+1)]<0): nn="nn+1" print("er="+str(E))
print(" nodes="+str(nn))
for i in range(m+1):
	print(" psir(%f)="%f"" % (x["x%s" i],p["p%s" i])) < code></0):>

这个代码实际上写的很拙劣，因为我目前只学了一点，各位不要嘲笑……

关于输入的一点说明：为了让得到的结果尽可能精确，$s$的选择应当尽量小一些，这样可以绘制更多的点。m为间隔的数量，显然$m=\frac{(x_{\text{max}}-x_0)}{s}$。

运行该代码，写入输入，则得到输出（同时会计算节点的数量）。

Other
开启自动换行复制
Enter initial xr:-4
Enter the increment sr:0.05
Enter the number of intervals m:160
Enter the reduced energy Er:0.5
Er=0.5
Nodes=0
Psir(-4.000000)=0.000000
Psir(-3.950000)=0.000100
Psir(-3.900000)=0.000204
Psir(-3.850000)=0.000315
...

其中最后一行为$x_{\text{max}}$的波函数值：

Other
开启自动换行复制
Psir(4.000000)=0.000489

可以看出，与0很接近。

若其偏离0，则需要手动对能量进行修改再次进行计算。

用以上数据在Origin中绘图：

与量子谐振子基态的波函数非常吻合（输出中表明了节点数为0）。

若修改初猜能量，得到第一激发态的波函数也是同样可行的。

但需要注意的一点是，当振动量子数为奇数时，$\psi$首先是小于零的。因此需要修改代码中的P["p1"]=0.0001为-0.0001。

P.S.书中提供了一段C++实现代码：
Numerov method.cpp 0.81KB CPP 81次下载

修改版代码

实际使用代码尝试后会发现$E_0,E_1,E_2,\dots$分别近似等于0.5，1.5，2.5...（不是严格相等大概是程序计算时数值精度的原因）

这是非常理所当然的一点，因为量子谐振子的能级满足：

$$E_v=(v+\frac{1}{2})h\nu$$

约化后能量自然等于0.5，1.5，2.5...

（以下内容我也不知道算不算画蛇添足）

所以我考虑通过修改代码，通过输入振动量子数，再利用程序自动对能量进行调整（由于程序原因，输入精确本征值，最后一点的波函数仍不够趋于0），最后进行数值计算。

所以我在进行最终计算前，设定一收敛条件，若$\psi (x_{\text{max}})$的绝对值小于0.0001，则判断$E_{\text{guess}}$足够准确。并由此进行计算。

修改后的代码如下（同求不嘲笑……）：

Other
开启自动换行复制
X={"x0":0.0}
P={"p0":0.0}
G={"g0":0.0}
nn=0

x=float(input("Enter initial xr:"))
s=float(input("Enter the increment sr:"))
m=int(input("Enter the number of intervals m:"))
v=float(input("Enter the vibration quantum number:"))

if v%2==0:
        P["p1"]=0.0001
else:
        P["p1"]=-0.0001

p1=P["p1"]
E=v+0.5
El=v
Eh=v+1
M=E
#定义函数求最后一点的psir
def iteration(E,x0,s,m,p1):
        p0=0.0
        x1=x0+s
        x2=x1+s
        ss=s*s/12
        for i in range(m-1):
                g0=x0*x0-2*E
                g1=x1+x1-2*E
                g2=x2*x2-2*E
                p2=(-p0+2*p1+10*g1*p1*ss+g0*p0*ss)/(1-g2*ss)
                x0=x0+s
                x1=x1+s
                x2=x2+s
                p0=p1
                p1=p2
                m=m-1
        return p2
#二分法
while abs(iteration(E,x,s,m,p1))>0.0001:
        if iteration(E,x,s,m,p1)>0:
                if iteration(Eh,x,s,m,p1)<0: e="(E+Eh)/2" if iteration(e,x,s,m,p1)>0:
                                pass
                        else:
                                Eh=E
                                E=M
                else:
                        E=(E+El)/2
                        if iteration(E,x,s,m,p1)>0:
                                pass
                        else:
                                El=E
                                E=M
        else:
                if iteration(Eh,x,s,m,p1)<0: e="(E+El)/2" if iteration(e,x,s,m,p1)<0: pass else: el="E" eh="E" x["x0"]="x" x["x1"]="X["x0"]+s" g["g0"]="X["x0"]*X["x0"]-2*E" g["g1"]="X["x1"]*X["x1"]-2*E" ss="s*s/12" for i in range(1,m): x["x%s" % (i+1)]="X["x%s"" i]+s g["g%s" (i+1)]*x["x%s" (i+1)]-2*e p["p%s" (i-1)]+2*p["p%s" i]+10*g["g%s" i]*p["p%s" i]*ss+g["g%s" (i-1)]*p["p%s" (i-1)]*ss) (1-g["g%s" (i+1)]*ss) (p["p%s" (i+1)]<0): nn="nn+1" print("er="+str(E))
print(" nodes="+str(nn))
for i in range(m+1):
        print(" psir(%f)="%f"" (x["x%s" i],p["p%s" i])) < code></0:></0:>

在上面的代码中，我定义了一个函数用于求$\psi (x_{\text{max}})$，再利用二分法调整能量$E$，当$\psi (x_{\text{max}})$满足收敛条件后，便使用调整后的能量计算各点波函数的数值（这部分与第一段代码没有区别）。

从原理上来说，这段代码应该没有明显的毛病。但是发现得到的能量并不准确。

测试后发现是因为函数的输出与第一段代码的输出有明显的差距。

把函数的代码单独运行，并使其输出每一次得到的$\psi$（前两行是我手动补上去的）

Other
开启自动换行复制
0.0000000000
0.0001000000
0.00019873416928030692
0.00029497854356514764
0.00038754990743566356
0.00047532063060121913
0.0005572325408706247
0.0006323096548955608
0.0006996696081080572
...
0.26353345349798035
0.3024049078455913
0.3474895062785628

这是第一段代码的输出：

Other
开启自动换行复制
Psir(-4.000000)=0.000000
Psir(-3.950000)=0.000100
Psir(-3.900000)=0.000204
Psir(-3.850000)=0.000315
Psir(-3.800000)=0.000436
Psir(-3.750000)=0.000573
Psir(-3.700000)=0.000728
...
Psir(3.900000)=0.000542
Psir(3.950000)=0.000506
Psir(4.000000)=0.000489

可以看出两段代码的输出结果的相差程度越来越大，到了最后一行，结果已经相差近三个数量级。用该函数的结果来调整能量的优化自然难以得到准确的结果。

我仔细检查了函数的代码，自认为没有错误。所以怀疑是浮点数精度（计算和储存时的四舍五入）的原因。

所以目前第二段代码暂不用来计算。