用松弛法解薛定諤方程

劉觀福，余 聰

(中山大學 物理與天文學院，廣東 珠海 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 更復雜勢能的情況

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