求解一維擴散方程的一種高精度緊致差分方法

楊曉佳，　葛永斌

(寧夏大學 數學計算機學院　寧夏 銀川 750021)

?

求解一維擴散方程的一種高精度緊致差分方法

楊曉佳，葛永斌

(寧夏大學 數學計算機學院寧夏 銀川 750021)

摘要：針對一維擴散方程，空間采用四階Padé 公式，時間采用廣義的梯形公式，差分離散得到了一種時間二階、空間四階精度的隱式緊致差分格式，其截斷誤差為O(τ2+h4). 通過理論分析證明了此格式是無條件穩定的. 最后通過數值實驗驗證了格式的精確性和可靠性.

關鍵詞：擴散方程； Padé 逼近； 緊致格式； 廣義梯形公式； 穩定性

0引言

擴散方程是一類重要的拋物型微分方程，它在數學、物理、化學等領域都有著廣泛的應用.自然環境、工程設備及生物機體中的許多物理現象都可用擴散方程來描述，由于物理問題本身的復雜性，其精確解往往不容易求得，因此研究其數值求解方法具有非常重要的意義.

求解該問題的方法已有很多，如有限差分法、有限元法和有限體積法等.在有限差分法研究方面，近年來，已有大量研究工作.如文獻[1]針對一維拋物型方程，采用中心差分格式將空間離散，得到一個關于時間t的半離散微分方程，然后采用廣義的梯形公式離散微分方程，得到了一個截斷誤差為O(τ2+h2)的隱式緊致差分格式.文獻[2]在離散關于時間t的半離散微分方程時采用了擴展的辛普森公式，得到了一個O(τ4+h2)的隱式緊致差分格式.文獻[3]采用含參數差分方程逼近的方法,構造了一維齊次擴散方程的高精度隱格式，隨參數選取的不同，格式的截斷誤差可分別達到O(τ2+h4)和O(τ4+h4)，但是O(τ4+h4)格式是無條件不穩定的. 文獻[4]基于泰勒級數展開和待定系數的方法,構造了一維拋物型方程精度為O(τ3+h4)的3層3點顯格式，其穩定性條件為r≤1/2(r為網格比)，即格式是條件穩定的.文獻[5]利用一階微商和二階微商的四階緊致差分逼近公式，推導出了求解一維擴散方程精度分別為O(τ2+h4)和O(τ4+h4)的兩種隱式差分格式，其中O(τ2+h4)格式是2層無條件穩定的，而O(τ4+h4)格式是3層無條件不穩定的. 文獻[6]提出了求解一維拋物型方程的一族2層6點隱式格式，格式的截斷誤差為O(τ2+h4)，并通過Fourier的方法證明了當1/2≤θ≤1(θ為參數)時，格式是無條件穩定的，當0≤θ<1/2，并且r≤1/6(1-2θ)，格式才是穩定的. 文獻[7]采用區域分裂的方法，構造了一種求解一維和二維拋物型方程的高精度顯格式，其截斷誤差為O(τ3+h3)，并且通過理論分析的方法證明其格式是無條件穩定的. 文獻[8]對空間變量應用中心差分格式和緊致差分格式離散，時間變量采用二級四階Runge-Kutta方法，構造了求解一維擴散方程的精度為O(τ4+h2)和O(τ4+h4)的兩種無條件穩定的隱式差分格式. 文獻[9]利用拋物型方程解的一階偏導數在特殊節點處的一個差分近似式和二階中心差商近似式，用待定系數法構造了一族隱式差分格式，格式的截斷誤差為O(τ3+h4)，穩定性條件為r>1/6. 文獻[10]針對一維拋物型方程的初邊值問題，采用待定系數法和泰勒級數展開的方法構造了一種2層8點隱式差分格式，其格式的截斷誤差為O(τ3+h5)，穩定性條件為0.001

1差分格式的建立

考慮一維含源項的擴散方程的初邊值問題:

ut=vuxx+f(x,t), 00,

(1)

給定初始條件為:

u(x,0)=φ(x), 0≤x≤l,

(2)

給定邊界條件為:

u(0,t)=a(t),u(l,t)=b(t),t>0,

(3)

其中：u(x,t)是待求未知函數；v(v>0)為擴散項系數；f(x,t)為非齊次項；a(t),b(t),φ(x)均為已知函數.

考慮在第n時刻的值，對空間內部節點采用四階精度的Padé公式[11]逼近，則有

(4)

對于空間邊界節點的處理，由式(1)與式(3)可得：

(5)

(6)

令

于是有

(7)

略去高階項，(7)式可變形為

(8)

其中：E=J-1S;C(tn)=J-1[C2(tn)-h2C1(tn)].

由(1)式和(8)式可得

(9)

對時間方向的離散，采用廣義梯形公式(GTF(a))[12]，有

(10)

(11)

將(9)式代入(10)式和(11)式可得：

(12)

(13)

將(12)式代入(13)式,整理可得

(14)

式(14)即為所構造的差分格式,通過對空間的離散過程不難發現，其具有四階精度. 另外，由文獻[1]的研究結論可知，時間具有二階精度，因此本文所構造的差分格式(14)的截斷誤差為O(τ2+h4).

2穩定性分析

下面分析上述差分格式(14)的穩定性，假定邊界條件精確滿足，且右端項f無誤差存在.

引理 1假設λ和x∈RN-1(x≠0)分別是矩陣J-1S的特征值和特征向量(J和S如上文定義)，則λ是實數且有λ>0.

證明由于λ和x分別是矩陣J-1S的特征值和特征向量，那么J-1Sx=λx?λxTJx=xTSx,令x=(x1,x2,…,xN-1)T,那么

所以xTSx>0， 而

則xTSx>0,由λxTJx=xTSx?λ>0，且λ為實數.

定理 1格式(14)是無條件穩定的.

證明令λj(j=1,2,…,N-1)是矩陣E的任意特征值，由引理1可知λj>0. 令

那么矩陣D的任意特征值為

記

則有

那么

因此，格式(14)是無條件穩定的.

3數值算例

為了驗證本文所提格式的精確性和可靠性，現考慮如下4個帶有精確解的問題，分別采用文獻[1,2，5，9]和本文格式(14)進行了計算，給出了不同時刻、不同時間步長、不同空間步長以及不同網格比r的最大絕對誤差及收斂階，其中最大絕對誤差及收斂階定義如下：

Error1和Error2分別為空間網格步長為h1和h2時的最大絕對誤差.

算例1[5]

ut=uxx, 00,

u(x,0)=sin(πx), 0≤x≤1,

u(0,t)=u(1,t)=0,t>0,

問題的精確解為u(x,t)=e-π2tsin(πx).

算例2[5]

ut=uxx+e-π2tsin(πx), 00,

u(x,0)=0, 0≤x≤1,

u(0,t)=u(1,t)=0,t>0,

問題的精確解為u(x,t)=te-π2tsin(πx).

算例3

算例4[1]

表1給出了問題1在t=1時對不同N的最大絕對誤差及收斂階，當問題1采用文獻[1]和本文格式(14)計算時，均取a=0.3.從表1可以看出，文獻[1]只有二階精度，并且其最大絕對誤差比本文格式的大幾個數量級，而文獻[5,9]和本文格式均具有四階精度，但本文格式的最大絕對誤差較文獻[5,9]的最大絕對誤差要小. 表2給出了N=32時，當網格比r分別為0.8、3.2、12.8、51.2時的最大絕對誤差，可以看出格式(14)是穩定的，這與本文的理論分析結果是一致的. 表3給出了問題2在t=0.25時，對不同N的最大絕對誤差及收斂階，用文獻[1]的方法和本文格式計算時，均取a=0.7.從表中的數值結果可以看出，本文格式比文獻[1]和文獻[5]的計算結果更為精確，并且當網格數逐漸增大時，本文格式的最大絕對誤差要比文獻[1]中的最大絕對誤差小幾個數量級. 表4給出了問題3的計算結果，該問題的邊界是關于時間t的函數，從表4中的計算結果可以看出，在推導本文格式的過程中對邊界條件的處理方法是切實可行的.對于該問題，用文獻[1]的方法和本文格式計算時，均取a=0.2. 表5給出了當N=32時，問題4在t=5時，文獻[1]、文獻[2]、文獻[5]和本文格式的數值解和精確解，在計算過程中，文獻[1]和本文格式中取a=1/3，同樣可以看出本文格式的計算結果較其他3個格式的計算結果更為精確. 表6給出了問題4當N=32時，當網格比r分別為0.8、3.2、12.8、51.2時的最大絕對誤差，從表中的數值結果同樣可知計算均是穩定的，這與本文的理論分析結果是一致的.

表1　問題1，當τ=h2,t=1時對不同N的最大絕對誤差及收斂階

表2　問題1，當N=32，t=0.5時對不同的r最大絕對誤差

表3　問題2，在t=0.25時對不同N的最大絕對誤差及收斂階

表4　問題3，在t=1時對不同N的最大絕對誤差及收斂階

表5　問題4，當N=32,τ=h2，t=5時的數值解和精確解

表6　問題4 當N=32，t=1時對不同r的最大絕對誤差

4結論

本文針對一維擴散方程提出了一種新的隱式高精度緊致差分格式，格式的截斷誤差為O(τ2+h4).通過理論分析證明了格式是無條件穩定的. 通過數值實驗將本文格式與其他文獻中的格式進行了比較，表明本文格式的計算結果更為精確，并且其精度和穩定性與理論分析結果相一致，從而驗證了本文格式的精確性與可靠性. 雖然本文方法只是針對一維問題的，但是本文的方法可以直接推廣到二維和三維，關于此方面的研究，我們將另文報道.

參考文獻：

[1]CHAWLA M M, Al-ZANAIDI M A, EVANS D J. Generalized trapezoidal formulas for parabolic equations[J]. International journal of computer mathematics,1999,70(30)：429—443.

[2]SALLAM S, NAIM A M, ABDEL-AZIZ M R. Unconditionally stable C1-cubic spline collocation method for solving parabolic equations[J]. International journal of computer mathematics,2004,81(7): 813—821.

[3]周順興. 解拋物型偏微分方程的高精度差分格式[J]. 計算數學,1982,4(2): 204—213.

[4]馬明書. 一維拋物型方程的一個新的高精度顯式差分格式[J]. 數值計算與計算機應用,2001,22(2): 156—160.

[5]葛永斌, 田振夫, 詹詠,等. 求解擴散方程的一種高精度隱式差分方法[J]. 上海理工大學學報,2005,27(2): 107—112.

[6]詹涌強. 求解拋物型方程的一族六點隱式差分格式[J]. 安徽大學學報(自然科學版),2012,36(4): 26—29.

[7]WEN R H, SHAO H Z. Domain decomposition schemes with high-order accuracy and unconditional stability[J]. Applied mathematics and computation,2013,219(11): 6170—6181.

[8]開依沙爾熱合曼, 努爾買買提黑力力. 求解擴散方程的二級四階隱式Runge-Kutta方法[J]. 湖北大學學報(自然科學版), 2014, 36(5): 476—480.

[9]詹涌強, 張傳林. 解拋物型方程的一族高精度隱式差分格式[J]. 應用數學和力學, 2014, 35(7): 790—797.

[10]周敏, 高學軍, 董超. 解拋物型方程的八點隱式差分格式[J]. 廣東工業大學學報(自然科學版), 2014, 31(4): 71—78.

[11]LELE S K. Compact finite difference schemes with spectral-like resolution[J]. Journal of computational physics,1992,103(1): 16—42.

[12]CHAWLA M M, Al-ZANAIDI M A, EVANS D J. A class of generalized trapezoidal formulas for the numercal integration of y′=f(x,y)[J]. International journal of computer mathematics, 1996, 62(1/2):131—142.

(責任編輯：方惠敏)

A High-order Compact Difference Method for the One-dimensional Diffusion Equation

YANG Xiaojia , GE Yongbin

(CollegeofMathematicsandComputerScience,NingxiaUniversity,Yinchuan750021,China)

Abstract:For the one dimensional diffusion equation, a high-order compact difference scheme was constructed by using the fourth-order Padé formula for space discretization and the generalized trapezoidal formula for time discretization, and the truncation error of the scheme was O(τ2+h4). The unconditional stability of the scheme was analyzed by theoritical method. Numerical experiments were carried out to verify the accuracy and the reliability of the present method.

Key words:diffusion equation; Padé approximation; compact scheme; generalized trapezoidal formula ; stability

收稿日期:2015-06-28

基金項目:國家自然科學基金資助項目(11361045, 11161036)；寧夏高等學校科學技術研究項目(NGY2013019).

作者簡介:楊曉佳(1988—)，女，寧夏吳忠人，碩士研究生，主要從事偏微分方程數值解法研究，E-mail：yang_xiaoj@sina.com；通信作者：葛永斌(1975—)，男，寧夏吳忠人，教授，博士生導師，主要從事偏微分方程數值解法研究，E-mail：gyb@nxu.edu.cn.

中圖分類號：O241.82

文獻標志碼：A

文章編號：1671-6841(2016)01-0010-07

DOI:10.3969/j.issn.1671-6841.201506013

引用本文：楊曉佳，葛永斌. 求解一維擴散方程的一種高精度緊致差分方法[J]. 鄭州大學學報(理學版)，2016，48(1)：10—16.