求解擴散方程的幾種差分格式的比較

祖麗胡瑪爾·卡迪爾，阿米娜·沙比爾

(喀什大學數學與統計學院，新疆 喀什 844000)

0 引言

擴散方程是一類基本的運動方程，是拋物型偏微分方程一個很重要的分支，在眾多領域都有著廣泛的應用，可用于環境科學、能源開發、流體力學和電子科學等許多領域[1].因此，數值求解拋物型偏微分方程問題具有重要的理論意義和實用價值.求解擴散方程的數值方法主要是有限差分法、有限元法、有限體積法等多種方法.但是對于對流占優問題，通常是用差分法或有限元法進行求解.本文利用三種差分格式求解了一維擴散方程的數值解，并給出了穩定性比較，最后通過數值算例進一步說明了數值解法的有效性.

1 差分格式的建立

考慮以下常系數擴散方程的初邊值問題：

其中，a>0為正常數，f(x,t),φ(x),α(t),β(t)都是已知函數，L為空間長度，T為時間長度.

首先，對求解區域D={(x,t)|0 ≤x≤L,0 ≤t≤T}作矩形網格剖分，將空間[0,L]等分m份，節點為xi=0+ih（0 ≤i≤m，h=）.再對時間區間[0,T]等分n份，節點為tj=0+jτ（0 ≤j≤n，.這樣得到(m+1)(n+1)個網格節點(xi,tj).

1.1 向前歐拉格式

在網格節點建立節點離散方程，即

關于時間的一階導數的向前差商形式為

上式誤差為O(τ).

類似地，關于空間的二階偏導數的中心差商形式為

上式誤差為O(h2)

將式（3）和（4）帶入（2）式中的第一式，再用數值解近似代替精確解u(xi,tj)，并忽略高階小項，即可得到向前歐拉格式：

經過整理可得向前歐拉格式：

由（6）式知，第j+1 個時間層上的數值解可由第j層上的已知信息表示出來.如此一層一層計算，可得到所有網格點的數值解，計算簡單、便捷[2].

1.2 向后歐拉格式

對時間的偏導數用向后差商，關于空間的二階偏導數用中心差商代替，再用數值解近似代替精確解u(xi,tj)，并忽略高階小項，即可得到向后歐拉格式：

實際計算時，在每個時間層都需要求解一個線性方程組，計算比較復雜[3].

1.3 Crank-Nicolson 差分格式

對時間的一階偏導數用中心差商代替，則有

上式誤差為O(τ2).

接下來處理虛擬節點處關于空間的二階偏導數

將上面兩式代入式（10）得

上式和（9）式一起代入（8）式，用數值解近似代替精確解u(xi,tj)，并忽略高階小項可得Crank-Nicolson差分格式：

易見，此格式的階段誤差為O(τ2+h2).

整理上述格式，可得

上式的矩陣形式可以利用追趕法來求解.

2 穩定性比較

向前歐拉格式[5]是條件穩定的，其穩定性條件為r=，如果網格比不滿足穩定性條件，導致最后結果失效.

向后歐拉格式是無條件穩定的，網格比很大情況下數值解仍然收斂.

Crank-Nicolson 差分格式[4]是二階無條件穩定的，其中一個至關重要的因素就是將原項方程弱化，使之在相鄰時間層網格節點的終點處成立，而不是在網格節點處成立.

3 數值算例

例：計算拋物型方程的初值問題

已知其精確解為u(x,t)=tsin(x).

由表1 和圖1 可見，r=0.5 時數值結果完全吻合；當網格步長不滿足穩定性條件時，誤差會以指數增長，結果不收斂，數值結果無效.這表明向前歐拉格式是條件穩定的.由表2 和圖2 可見，無論r的取值如何，也就是無論時間、空間步長如何選取，向后歐拉格式恒穩定，即該格式是無條件穩定的.由表3 和圖3 可見，Crank-Nicolson 格式是無條件穩定的，差分格式是二階精度的.

表1 當r=0.5時，向前歐拉格式的數值解和精確解

表2 向后歐拉格式在不同網格步長下的數值解

表3 Crank-Nicolson差分格式在不同網格步長下的數值解

圖1 向前歐拉格式

圖2 向后前歐拉格式

圖3 Crank-Nicolson格式

4 結論

本文討論了求解拋物型方程初邊值問題的幾種常用方法，如向前歐拉格式方法、向后歐拉格式方法、Crank-Nicolson 格式方法，并了解各種方法的優缺點.結果表明，向前歐拉格式方法計算簡單、直接，但穩定性差，只有在時間、空間步長的合理選取之下才能保證數值解的收斂性；向后歐拉格式方法則彌補了向前歐拉格式方法的缺點，不需要考慮時間、空間步長的選取，但每個時間層上都要求解一個線性方程組，計算量會加大.這兩個格式方法的精度都是時間一階、空間二階；Crank-Nicolson 格式方法是無條件穩定的，時間和空間都是二階精度的.