劉觀福，余 聰

(中山大學 物理與天文學院，廣東 珠海 519082)

在經(jīng)典情況下，最重要的常微分方程是牛頓第二定律的方程.在量子情況下，最重要的方程的是薛定諤方程.除了一些比較的簡單的情況，這兩個常微分方程通常難以找到解析解的，必須用數(shù)值解法，通過數(shù)值解來分析.一般情況下，牛頓第二定律的方程和薛定諤方程在數(shù)值解法上是有區(qū)別的.從已知的條件來說，對于經(jīng)典情況的物體，它的初始位置、初始速度以及各處的加速度是明確的，解它對應的牛頓第二定律方程是一個初值問題；而對量子情況下的粒子，一般只明確它在邊界上的概率幅函數(shù)的值，對應的薛定諤方程是一個邊值問題.從最終數(shù)值解法的結果來說，牛頓第二定律方程一般只有一個解，薛定諤方程往往有無數(shù)個解.

1 打靶法和松弛法的區(qū)別

常微分方程的數(shù)值解法的思想是將微分方程改寫成有限差分方程[1-8].打靶法和松弛法都是采用了這一思想.它們的區(qū)別在于如何去處理邊界條件.物理中很常見的情況是知道初始點和終點的函數(shù)值，對于一維無限深勢阱情況下的概率幅函數(shù)，在兩個邊界上的概率幅函數(shù)值是0[2].打靶法的基礎是龍格庫塔方法，需要得到一個邊界點上的函數(shù)值和函數(shù)的導數(shù)值.在導數(shù)值不知道的情況下，采取的解決辦法是在一定范圍內猜測導數(shù)值，從初始位置邊界點開始計算，計算到終止位置邊界點，然后比較終止位置邊界點計算的函數(shù)值和預期值，如果計算值和預期值在誤差允許范圍內，就結束計算，反之，就重新猜測初始位置的導數(shù)值，重新計算，直到符合誤差允許范圍.這個過程和打靶射擊是極其相似的，所以也被稱為“打靶法”[1,8].而松弛法是一開始猜測在要計算的區(qū)間內猜測函數(shù)值和導數(shù)值，然后根據(jù)差分方程和邊界值來計算偏差值，偏差值符合要求就結束計算，反之減去偏差值，再次計算偏差值，直到迭代出符合要求的解.這個過程就像一個繃緊的橡皮筋，一步步地松弛到最理想的狀態(tài)，所以稱為“松弛法”.打靶法是把一個邊界當成計算結果是否符合要求的判據(jù)，在一次計算過程中，這個邊界是不起作用的，只在一次計算完成后用來判斷這次計算是否符合要求；松弛法是在每次迭代中都要用到所有邊界條件來計算偏差.從這層面上來理解，松弛法更有效地利用了邊界條件，計算速度會更快.實際上，打靶法猜測初始的導數(shù)值會在根本上影響一次計算出來的結果，就像打靶的時候，初始的射擊角度決定了最終擊中的位置，所以需要反復多次的去調整，才能得到符合要求的結果，所以需要的計算時間長；而松弛法，雖然也要猜測初始值，但是初始值猜測不夠準確，還可以通過迭代修正，就像一個橡皮筋，初始繃得很緊，但是可以在后面一步步松弛到理想狀態(tài)，所以需要的計算時間短.

2 松弛法的計算原理

2.1 基本算法

考慮N個相互耦合的一階常微分方程，一個N階常微分方程可以轉換為N個一階常微分方程.對于一個一般的一階常微分方程：

(1)

將要計算的區(qū)間均勻分成M-1份，得到包括邊界點在內的M個格點，每個格點用k來區(qū)分，然后將常微分方程改寫成有限差分方程：

(2)

考慮到有N個待求的未知函數(shù)，定義有限差分方程的誤差函數(shù)為

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

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

(3)

上面由于N個未知函數(shù)耦合在一起，所以式(2)中的y要改成N維向量.注意到上面的誤差函數(shù)只是中間點的誤差函數(shù)，這表示中間點和有限差分方程的偏差.下面定義邊界點的誤差函數(shù).

設初始位置有n1個邊界條件，則初始位置邊界的誤差函數(shù)：

E1=B(x1,y1)

(4)

終止位置有N-n1個邊界條件，則終止位置邊界的誤差函數(shù)：

EM+1=C(x1,yM)

(5)

這里N個邊界條件不集中在一個邊界上，因為如果集中在一個邊界上，那么直接用龍格庫塔方法就可以求解，不需要打靶法和松弛法.

到這里，我們可以由式(3)～(5)計算出每一次迭代的N個函數(shù)在M個點的函數(shù)值和邊界條件以及有限差分方程的誤差.接下來需要建立誤差函數(shù)和每一次迭代的函數(shù)修正量Δy直接的關系.公式里面

(6)

所以每一個xk是確定的，誤差函數(shù)Ek只是yk的函數(shù).修正之后的誤差函數(shù)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(tǒng)+Δy，然后再將誤差函數(shù)和設定的精度比較，不符合要求，就再次計算修正量Δy，繼續(xù)迭代，直到符合精度要求.

2.2 整體的相對誤差、誤差范數(shù)、校正參數(shù)

在1.1中介紹了基本迭代方法.下面介紹誤差控制的方法.定義整體的相對誤差err

(10)

式(10)的scalv(j)，(j=1,2,...N)稱之為誤差范數(shù)[1,8].誤差范數(shù)是事先給定的，用來大致描述N個不同的待求解的函數(shù)的量級.因為N個不同的待求解函數(shù)有各自量級，對應的|Δyj,k|量級也不一樣，單純把它們加起來，會導致量級大的函數(shù)占主導，而量級小的函數(shù)的收斂性差.除以各自的誤差范數(shù)相當于把各個待求解的函數(shù)的相對誤差加起來了，然后再除以個數(shù)，得到相對誤差的平均值err，稱之為整體的相對誤差.計算的時候期望err符合精度要求，符合精度要求就停止計算，返回函數(shù)值.

在迭代的過程中，1.1中提及的將y+Δy去代替y，這里引入校正參數(shù)slowc[1,8]：

(11)

當整體的相對誤差err比較大的時候，迭代的時候就用修正量Δy的一部分加上y，得到迭代后的y，整體的相對誤差小的時候，就用修正量Δy加上y，得到迭代后的y.這個也不難理解：向上緊繃的橡皮筋在松弛的過程中可能瞬間變成向下的緊繃狀態(tài)，在松弛的過程中出現(xiàn)了修正過量的情況，出現(xiàn)這種情況，我們需要削減修正量.

3 無限深勢阱下的薛定諤方程

3.1 薛定諤方程

無限深勢阱的邊界如圖1所示.

圖1 無限深勢阱的勢能圖

對應的薛定諤方程

(12)

勢能函數(shù)

(13)

所以ψ(x)滿足

(14)

所以計算的區(qū)間是[-1,1]，邊界條件是-1和+1處ψ(x)=0.除此之外，還要滿足歸一化條件

(15)

3.2 誤差函數(shù)

為了方便敘述，引入3個函數(shù)：

(16)

這3個函數(shù)滿足

(17)

邊界條件

(18)

這里面把ψ(x)在+1處的導數(shù)值設定為5，這里是為了計算的方便.式(15)不方便寫出誤差函數(shù).同時因為式(15)是用來限制ψ(x)的范圍，這個條件改成ψ′(1)=5，在計算結果上只是把ψ(x)變成了cψ(x)(c是一個常數(shù)).只需要最后再把計算出來的cψ(x)乘以c的倒數(shù)就得到了ψ(x).在后文的計算結果的展示中，未做乘以c的倒數(shù)這一步操作.和前文保持一致，將區(qū)間[-1,1]分成包括邊界在內的M個計算格點.

中間M-2個格點的誤差函數(shù)：

(19)

對應的矩陣S的元素：

(20)

(21)

(22)

邊界點x=-1處的誤差函數(shù)：

E3,1=y1,1

(23)

對應的矩陣S的元素：

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

(24)

邊界點x=1處的誤差函數(shù)：

(25)

對應的矩陣S的元素：

(26)

3.3 理論分析

由式(16)～(18),不難得到未歸一化后的理論解(沒有乘以c的倒數(shù)的理論解，概率幅平方積分不是1)

(27)

3.4 理論和計算的對比

就像橡皮筋松弛之后的狀態(tài)依賴于初始緊繃的狀態(tài)一樣，松弛法計算的結果依賴于初始值的設定.一般來說，松弛法計算的結果和初始值之間的關系是不明顯的.但是對應無限深勢阱的薛定諤方程，還是有一定的規(guī)律可循.在不做任何理論解分析的情況下，可以先隨機猜測一個邊界處ψ(x)的導數(shù)值，然后用龍格-庫塔法去計算[1]，雖然這種計算往往不能滿足另一個邊界的條件，但是從ψ(x)的圖像上可以發(fā)現(xiàn)，ψ(x)具有很好的周期性，甚至可以進一步猜測是正弦函數(shù).這里，猜測ψ(x)和ψ′(x)的初始值是方波型函數(shù).k2的初始猜測值是一個正數(shù)即可.方波函數(shù)的周期滿足

(28)

下面是計算情況的展示，如圖2～10所示.計算時設定的整體的相對誤差err是5×10-5，也就是說最終計算結果滿足err≤5×10-5.計算的時候用的是Fortran語言，畫圖用的是Python語言，圖上Fortran語言的線是數(shù)值計算的結果.圖4中k2的理論解和數(shù)值解十分接近，圖上的縱坐標的單位長度是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 無限深勢阱的勢能圖

勢能函數(shù)：

(29)

一般會做如下代換[4]：

(30)

所以ψ(ξ)滿足

(31)

計算機計算的區(qū)間不能是[-∞,+∞]這里面計算的區(qū)間取為[-10,10]，即邊界條件是-10和+10處ψ(ξ)=0.同前文，將歸一化條件改為

ψ′(ξ)=-24

(32)

4.2 理論和計算的對比

結果如圖12—圖14所示.

圖12 ψ(x)的初始猜測值和計算結果對比

圖13 ψ′(x)的初始猜測值和計算結果對比

圖14 λ-ξ2的初始猜測值和計算結果對比

5 計算結果評估

5.1 初始猜測值和最終結果之間的關系

從前面的圖可以看出：在無限深勢阱的情況下：初始猜測值和最后的數(shù)值解具有相同的周期性.猜測的方波函數(shù)和計算出來的圖像的周期性一致.這可以從傅立葉變換的角度來理解，方波函數(shù)傅立葉變換后，主要成分是和它周期一致的正弦函數(shù)，最后計算也收斂到這個正弦函數(shù).

5.2 計算結果的精確度

從圖12—圖14上看，計算結果和理論結果符合得較好，沒有太大的誤差.計算使用的程序設定的整體相對誤差要求是err≤5×10-5.這也說明計算的結果誤差較小.

5.3 計算耗時

區(qū)間分為4500個格點，一次編譯加運行的耗時需要5秒左右，耗時較少.

6 和其它方法的對比

如表1所示.

表1 各種情況下松弛法和其他方法的對比

6.1 無限或有限深勢阱的情況

無限或有限深勢阱的情況比較簡單.無限深勢阱下的薛定諤方程最后歸結為解二階常系數(shù)齊次線性微分方程.二階常系數(shù)齊次線性微分方程的兩種基礎解系是指數(shù)型函數(shù)(正弦和余弦函數(shù)可以用指數(shù)函數(shù)表示)[3].解出這個方程的一般形式的解，然后再令這個解符合邊界條件，就解出了這種情況下的薛定諤方程[2].

和松弛法對比起來，直接去解這個微分方程的優(yōu)勢在于不需要太大的計算量，同時解完后這個微分方程的所有解都出來了.對于松弛法來說，通過簡單地修改初始猜測值，通過計算機計算，可以很快拿到幾個該微分方程的解，然后通過畫圖，可以很快地對這個方程的解畫圖，得到解的性質.無論是從高校教師教學的角度，還是大學生學習的角度，將這兩種方法結合起來可以加深對薛定諤方程的理解，也可以更好地學習數(shù)值方法.

6.2 諧振子勢的情況

諧振子勢的情況比無限深勢阱的情況更加復雜.諧振子勢的情況下的薛定諤方程是二階變系數(shù)線性微分方程，無限深勢阱的是二階常系數(shù)線性微分方程[3].解諧振子勢的薛定諤方程方法，筆者目前知道的有：一種是先分析漸近行為，得到漸近行為情況下的解，然后一般情況的解先猜測是漸近解乘以一個多項式，代回薛定諤方程里面，比較兩邊的系數(shù)，得到多項式的具體形式[4]；另一種是通過構造哈密頓量與諧振子系統(tǒng)哈密頓量對易的超對稱系統(tǒng)，量子諧振子的性質就可以通過對超對稱系統(tǒng)的研究來得到[5,7].這兩種方法本質上都是利用了級數(shù)展開的思想，通過繁瑣代數(shù)運算得到解.

從無限深勢阱到諧振子勢，在代數(shù)上是兩種不同的問題，在解法上也很不一樣，運算的繁瑣程度大幅度提高了.而從無限深勢阱到諧振子勢，對于松弛法來說，只需要稍微改一下計算無限深勢阱時輸入的系數(shù)矩陣參數(shù)，就可以計算諧振子勢的薛定諤方程的解.這個角度來說，松弛法對不同勢能情況下的薛定諤方程的適用性更強.

6.3 更復雜勢能的情況

在勢能函數(shù)更加復雜的情況下，代數(shù)方法已經(jīng)很難得到準確解.即使得到了準確解，準確解的形式往往很復雜，不容易去分析準確解的增減性、周期性等等的性質.而實際應用中，并不是一定需要得到準確解，往往只需要得到一個符合要求的精確解.對符合要求的精確解進行分析就可以滿足需要.這時候，松弛法就顯得尤為重要.因為解簡單勢能問題的松弛法的代碼，稍微加以修改初始猜測值、計算區(qū)間、誤差函數(shù)等等參數(shù)，就可以用來解更復雜勢能下的問題.所以松弛法具有比代數(shù)方法更強的適用性.