用松弛法解薛定谔方程

刘观福，余 聪

(中山大学 物理与天文学院，广东 珠海 519082)

在经典情况下，最重要的常微分方程是牛顿第二定律的方程.在量子情况下，最重要的方程的是薛定谔方程.除了一些比较的简单的情况，这两个常微分方程通常难以找到解析解的，必须用数值解法，通过数值解来分析.一般情况下，牛顿第二定律的方程和薛定谔方程在数值解法上是有区别的.从已知的条件来说，对于经典情况的物体，它的初始位置、初始速度以及各处的加速度是明确的，解它对应的牛顿第二定律方程是一个初值问题；而对量子情况下的粒子，一般只明确它在边界上的概率幅函数的值，对应的薛定谔方程是一个边值问题.从最终数值解法的结果来说，牛顿第二定律方程一般只有一个解，薛定谔方程往往有无数个解.

1 打靶法和松弛法的区别

常微分方程的数值解法的思想是将微分方程改写成有限差分方程[1-8].打靶法和松弛法都是采用了这一思想.它们的区别在于如何去处理边界条件.物理中很常见的情况是知道初始点和终点的函数值，对于一维无限深势阱情况下的概率幅函数，在两个边界上的概率幅函数值是0[2].打靶法的基础是龙格库塔方法，需要得到一个边界点上的函数值和函数的导数值.在导数值不知道的情况下，采取的解决办法是在一定范围内猜测导数值，从初始位置边界点开始计算，计算到终止位置边界点，然后比较终止位置边界点计算的函数值和预期值，如果计算值和预期值在误差允许范围内，就结束计算，反之，就重新猜测初始位置的导数值，重新计算，直到符合误差允许范围.这个过程和打靶射击是极其相似的，所以也被称为“打靶法”[1,8].而松弛法是一开始猜测在要计算的区间内猜测函数值和导数值，然后根据差分方程和边界值来计算偏差值，偏差值符合要求就结束计算，反之减去偏差值，再次计算偏差值，直到迭代出符合要求的解.这个过程就像一个绷紧的橡皮筋，一步步地松弛到最理想的状态，所以称为“松弛法”.打靶法是把一个边界当成计算结果是否符合要求的判据，在一次计算过程中，这个边界是不起作用的，只在一次计算完成后用来判断这次计算是否符合要求；松弛法是在每次迭代中都要用到所有边界条件来计算偏差.从这层面上来理解，松弛法更有效地利用了边界条件，计算速度会更快.实际上，打靶法猜测初始的导数值会在根本上影响一次计算出来的结果，就像打靶的时候，初始的射击角度决定了最终击中的位置，所以需要反复多次的去调整，才能得到符合要求的结果，所以需要的计算时间长；而松弛法，虽然也要猜测初始值，但是初始值猜测不够准确，还可以通过迭代修正，就像一个橡皮筋，初始绷得很紧，但是可以在后面一步步松弛到理想状态，所以需要的计算时间短.

2 松弛法的计算原理

2.1 基本算法

考虑N个相互耦合的一阶常微分方程，一个N阶常微分方程可以转换为N个一阶常微分方程.对于一个一般的一阶常微分方程：

(1)

将要计算的区间均匀分成M-1份，得到包括边界点在内的M个格点，每个格点用k来区分，然后将常微分方程改写成有限差分方程：

(2)

考虑到有N个待求的未知函数，定义有限差分方程的误差函数为

Ek≡yk-yk-1-(xk-xk-1)g(xk,xk-1,yk,yk-1),

(k=2,3,...,M)

(3)

上面由于N个未知函数耦合在一起，所以式(2)中的y要改成N维向量.注意到上面的误差函数只是中间点的误差函数，这表示中间点和有限差分方程的偏差.下面定义边界点的误差函数.

设初始位置有n1个边界条件，则初始位置边界的误差函数：

E1=B(x1,y1)

(4)

终止位置有N-n1个边界条件，则终止位置边界的误差函数：

EM+1=C(x1,yM)

(5)

这里N个边界条件不集中在一个边界上，因为如果集中在一个边界上，那么直接用龙格库塔方法就可以求解，不需要打靶法和松弛法.

到这里，我们可以由式(3)～(5)计算出每一次迭代的N个函数在M个点的函数值和边界条件以及有限差分方程的误差.接下来需要建立误差函数和每一次迭代的函数修正量Δy直接的关系.公式里面

(6)

所以每一个xk是确定的，误差函数Ek只是yk的函数.修正之后的误差函数Ek：

Ek(yk+Δyk,yk-1+Δyk-1)≈Ek(yk,yk-1)+

(k=1,3,...,M+1)

(7)

我们希望修正之后的Ek(yk+Δyk,yk-1+Δyk-1)等于0，所以得到

(k=1,2,...,M+1)

(8)

把式(8)改写成矩阵形式，约等号改为等于号,则

(j=1,2,...,N)，(k=1,2,...,M+1),

(9)

每一次迭代都用式(9)计算出修正量Δy，得到y+Δy，然后再将误差函数和设定的精度比较，不符合要求，就再次计算修正量Δy，继续迭代，直到符合精度要求.

2.2 整体的相对误差、误差范数、校正参数

在1.1中介绍了基本迭代方法.下面介绍误差控制的方法.定义整体的相对误差err

(10)

式(10)的scalv(j)，(j=1,2,...N)称之为误差范数[1,8].误差范数是事先给定的，用来大致描述N个不同的待求解的函数的量级.因为N个不同的待求解函数有各自量级，对应的|Δyj,k|量级也不一样，单纯把它们加起来，会导致量级大的函数占主导，而量级小的函数的收敛性差.除以各自的误差范数相当于把各个待求解的函数的相对误差加起来了，然后再除以个数，得到相对误差的平均值err，称之为整体的相对误差.计算的时候期望err符合精度要求，符合精度要求就停止计算，返回函数值.

在迭代的过程中，1.1中提及的将y+Δy去代替y，这里引入校正参数slowc[1,8]：

(11)

当整体的相对误差err比较大的时候，迭代的时候就用修正量Δy的一部分加上y，得到迭代后的y，整体的相对误差小的时候，就用修正量Δy加上y，得到迭代后的y.这个也不难理解：向上紧绷的橡皮筋在松弛的过程中可能瞬间变成向下的紧绷状态，在松弛的过程中出现了修正过量的情况，出现这种情况，我们需要削减修正量.

3 无限深势阱下的薛定谔方程

3.1 薛定谔方程

无限深势阱的边界如图1所示.

图1 无限深势阱的势能图

对应的薛定谔方程

(12)

势能函数

(13)

所以ψ(x)满足

(14)

所以计算的区间是[-1,1]，边界条件是-1和+1处ψ(x)=0.除此之外，还要满足归一化条件

(15)

3.2 误差函数

为了方便叙述，引入3个函数：

(16)

这3个函数满足

(17)

边界条件

(18)

这里面把ψ(x)在+1处的导数值设定为5，这里是为了计算的方便.式(15)不方便写出误差函数.同时因为式(15)是用来限制ψ(x)的范围，这个条件改成ψ′(1)=5，在计算结果上只是把ψ(x)变成了cψ(x)(c是一个常数).只需要最后再把计算出来的cψ(x)乘以c的倒数就得到了ψ(x).在后文的计算结果的展示中，未做乘以c的倒数这一步操作.和前文保持一致，将区间[-1,1]分成包括边界在内的M个计算格点.

中间M-2个格点的误差函数：

(19)

对应的矩阵S的元素：

(20)

(21)

(22)

边界点x=-1处的误差函数：

E3,1=y1,1

(23)

对应的矩阵S的元素：

S3,4=1,S3,5=0,S3,6=0

(24)

边界点x=1处的误差函数：

(25)

对应的矩阵S的元素：

(26)

3.3 理论分析

由式(16)～(18),不难得到未归一化后的理论解(没有乘以c的倒数的理论解，概率幅平方积分不是1)

(27)

3.4 理论和计算的对比

就像橡皮筋松弛之后的状态依赖于初始紧绷的状态一样，松弛法计算的结果依赖于初始值的设定.一般来说，松弛法计算的结果和初始值之间的关系是不明显的.但是对应无限深势阱的薛定谔方程，还是有一定的规律可循.在不做任何理论解分析的情况下，可以先随机猜测一个边界处ψ(x)的导数值，然后用龙格-库塔法去计算[1]，虽然这种计算往往不能满足另一个边界的条件，但是从ψ(x)的图像上可以发现，ψ(x)具有很好的周期性，甚至可以进一步猜测是正弦函数.这里，猜测ψ(x)和ψ′(x)的初始值是方波型函数.k2的初始猜测值是一个正数即可.方波函数的周期满足

(28)

下面是计算情况的展示，如图2～10所示.计算时设定的整体的相对误差err是5×10-5，也就是说最终计算结果满足err≤5×10-5.计算的时候用的是Fortran语言，画图用的是Python语言，图上Fortran语言的线是数值计算的结果.图4中k2的理论解和数值解十分接近，图上的纵坐标的单位长度是10-8.图片上看上去符合得不好，实际上精度很高，后面的图7和10也是同理.

图2 ψ(x)的初始猜测值和计算结果对比

图3 ψ′(x)的初始猜测值和计算结果对比

图4 k2的初始猜测值和计算结果对比

图5 ψ(x)的初始猜测值和计算结果对比

图6 ψ′(x)的初始猜测值和计算结果对比

当n=1,T=4时：

图7 k2的初始猜测值和计算结果对比

当n=2,T=2时：

图8 ψ(x)的初始猜测值和计算结果对比

图9 ψ′(x)的初始猜测值和计算结果对比

图10 k2的初始猜测值和计算结果对比

4 谐振子势下的薛定谔方程

4.1 薛定谔方程

谐振子势如图11所示.

图11 无限深势阱的势能图

势能函数：

(29)

一般会做如下代换[4]：

(30)

所以ψ(ξ)满足

(31)

计算机计算的区间不能是[-∞,+∞]这里面计算的区间取为[-10,10]，即边界条件是-10和+10处ψ(ξ)=0.同前文，将归一化条件改为

ψ′(ξ)=-24

(32)

4.2 理论和计算的对比

结果如图12—图14所示.

图12 ψ(x)的初始猜测值和计算结果对比

图13 ψ′(x)的初始猜测值和计算结果对比

图14 λ-ξ2的初始猜测值和计算结果对比

5 计算结果评估

5.1 初始猜测值和最终结果之间的关系

从前面的图可以看出：在无限深势阱的情况下：初始猜测值和最后的数值解具有相同的周期性.猜测的方波函数和计算出来的图像的周期性一致.这可以从傅立叶变换的角度来理解，方波函数傅立叶变换后，主要成分是和它周期一致的正弦函数，最后计算也收敛到这个正弦函数.

5.2 计算结果的精确度

从图12—图14上看，计算结果和理论结果符合得较好，没有太大的误差.计算使用的程序设定的整体相对误差要求是err≤5×10-5.这也说明计算的结果误差较小.

5.3 计算耗时

区间分为4500个格点，一次编译加运行的耗时需要5秒左右，耗时较少.

6 和其它方法的对比

如表1所示.

表1 各种情况下松弛法和其他方法的对比

6.1 无限或有限深势阱的情况

无限或有限深势阱的情况比较简单.无限深势阱下的薛定谔方程最后归结为解二阶常系数齐次线性微分方程.二阶常系数齐次线性微分方程的两种基础解系是指数型函数(正弦和余弦函数可以用指数函数表示)[3].解出这个方程的一般形式的解，然后再令这个解符合边界条件，就解出了这种情况下的薛定谔方程[2].

和松弛法对比起来，直接去解这个微分方程的优势在于不需要太大的计算量，同时解完后这个微分方程的所有解都出来了.对于松弛法来说，通过简单地修改初始猜测值，通过计算机计算，可以很快拿到几个该微分方程的解，然后通过画图，可以很快地对这个方程的解画图，得到解的性质.无论是从高校教师教学的角度，还是大学生学习的角度，将这两种方法结合起来可以加深对薛定谔方程的理解，也可以更好地学习数值方法.

6.2 谐振子势的情况

谐振子势的情况比无限深势阱的情况更加复杂.谐振子势的情况下的薛定谔方程是二阶变系数线性微分方程，无限深势阱的是二阶常系数线性微分方程[3].解谐振子势的薛定谔方程方法，笔者目前知道的有：一种是先分析渐近行为，得到渐近行为情况下的解，然后一般情况的解先猜测是渐近解乘以一个多项式，代回薛定谔方程里面，比较两边的系数，得到多项式的具体形式[4]；另一种是通过构造哈密顿量与谐振子系统哈密顿量对易的超对称系统，量子谐振子的性质就可以通过对超对称系统的研究来得到[5,7].这两种方法本质上都是利用了级数展开的思想，通过繁琐代数运算得到解.

从无限深势阱到谐振子势，在代数上是两种不同的问题，在解法上也很不一样，运算的繁琐程度大幅度提高了.而从无限深势阱到谐振子势，对于松弛法来说，只需要稍微改一下计算无限深势阱时输入的系数矩阵参数，就可以计算谐振子势的薛定谔方程的解.这个角度来说，松弛法对不同势能情况下的薛定谔方程的适用性更强.

6.3 更复杂势能的情况

在势能函数更加复杂的情况下，代数方法已经很难得到准确解.即使得到了准确解，准确解的形式往往很复杂，不容易去分析准确解的增减性、周期性等等的性质.而实际应用中，并不是一定需要得到准确解，往往只需要得到一个符合要求的精确解.对符合要求的精确解进行分析就可以满足需要.这时候，松弛法就显得尤为重要.因为解简单势能问题的松弛法的代码，稍微加以修改初始猜测值、计算区间、误差函数等等参数，就可以用来解更复杂势能下的问题.所以松弛法具有比代数方法更强的适用性.