劉觀福，余 聰

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

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

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

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

2 松弛法的計(jì)算原理

2.1 基本算法

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

(1)

將要計(jì)算的區(qū)間均勻分成M-1份，得到包括邊界點(diǎn)在內(nèi)的M個格點(diǎn)，每個格點(diǎn)用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ù)只是中間點(diǎn)的誤差函數(shù)，這表示中間點(diǎn)和有限差分方程的偏差.下面定義邊界點(diǎn)的誤差函數(shù).

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

E1=B(x1,y1)

(4)

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

EM+1=C(x1,yM)

(5)

這里N個邊界條件不集中在一個邊界上，因?yàn)槿绻性谝粋€邊界上，那么直接用龍格庫塔方法就可以求解，不需要打靶法和松弛法.

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

(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)計(jì)算出修正量Δy，得到y(tǒng)+Δy，然后再將誤差函數(shù)和設(shè)定的精度比較，不符合要求，就再次計(jì)算修正量Δy，繼續(xù)迭代，直到符合精度要求.

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

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

(10)

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

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

(11)

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

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

3.1 薛定諤方程

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

圖1 無限深勢阱的勢能圖

對應(yīng)的薛定諤方程

(12)

勢能函數(shù)

(13)

所以ψ(x)滿足

(14)

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

(15)

3.2 誤差函數(shù)

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

(16)

這3個函數(shù)滿足

(17)

邊界條件

(18)

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

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

(19)

對應(yīng)的矩陣S的元素：

(20)

(21)

(22)

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

E3,1=y1,1

(23)

對應(yīng)的矩陣S的元素：

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

(24)

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

(25)

對應(yīng)的矩陣S的元素：

(26)

3.3 理論分析

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

(27)

3.4 理論和計(jì)算的對比

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

(28)

下面是計(jì)算情況的展示，如圖2～10所示.計(jì)算時設(shè)定的整體的相對誤差err是5×10-5，也就是說最終計(jì)算結(jié)果滿足err≤5×10-5.計(jì)算的時候用的是Fortran語言，畫圖用的是Python語言，圖上Fortran語言的線是數(shù)值計(jì)算的結(jié)果.圖4中k2的理論解和數(shù)值解十分接近，圖上的縱坐標(biāo)的單位長度是10-8.圖片上看上去符合得不好，實(shí)際上精度很高，后面的圖7和10也是同理.

圖2 ψ(x)的初始猜測值和計(jì)算結(jié)果對比

圖3 ψ′(x)的初始猜測值和計(jì)算結(jié)果對比

圖4 k2的初始猜測值和計(jì)算結(jié)果對比

圖5 ψ(x)的初始猜測值和計(jì)算結(jié)果對比

圖6 ψ′(x)的初始猜測值和計(jì)算結(jié)果對比

當(dāng)n=1,T=4時：

圖7 k2的初始猜測值和計(jì)算結(jié)果對比

當(dāng)n=2,T=2時：

圖8 ψ(x)的初始猜測值和計(jì)算結(jié)果對比

圖9 ψ′(x)的初始猜測值和計(jì)算結(jié)果對比

圖10 k2的初始猜測值和計(jì)算結(jié)果對比

4 諧振子勢下的薛定諤方程

4.1 薛定諤方程

諧振子勢如圖11所示.

圖11 無限深勢阱的勢能圖

勢能函數(shù)：

(29)

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

(30)

所以ψ(ξ)滿足

(31)

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

ψ′(ξ)=-24

(32)

4.2 理論和計(jì)算的對比

結(jié)果如圖12—圖14所示.

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

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

圖14 λ-ξ2的初始猜測值和計(jì)算結(jié)果對比

5 計(jì)算結(jié)果評估

5.1 初始猜測值和最終結(jié)果之間的關(guān)系

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

5.2 計(jì)算結(jié)果的精確度

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

5.3 計(jì)算耗時

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

6 和其它方法的對比

如表1所示.

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

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

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

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

6.2 諧振子勢的情況

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

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

6.3 更復(fù)雜勢能的情況

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