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

楊曉佳，　葛永斌

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

?

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

楊曉佳，葛永斌

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

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

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

0引言

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

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

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

(4)

對于空間邊界節(jié)點的處理，由式(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)即為所構造的差分格式,通過對空間的離散過程不難發(fā)現(xiàn)，其具有四階精度. 另外，由文獻[1]的研究結論可知，時間具有二階精度，因此本文所構造的差分格式(14)的截斷誤差為O(τ2+h4).

2穩(wěn)定性分析

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

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

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

所以xTSx>0， 而

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

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

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

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

記

則有

那么

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

3數(shù)值算例

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

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

表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時的數(shù)值解和精確解

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

4結論

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

參考文獻：

[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]. 計算數(shù)學,1982,4(2): 204—213.

[4]馬明書. 一維拋物型方程的一個新的高精度顯式差分格式[J]. 數(shù)值計算與計算機應用,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]. 應用數(shù)學和力學, 2014, 35(7): 790—797.

[10]周敏, 高學軍, 董超. 解拋物型方程的八點隱式差分格式[J]. 廣東工業(yè)大學學報(自然科學版), 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—)，女，寧夏吳忠人，碩士研究生，主要從事偏微分方程數(shù)值解法研究，E-mail：yang_xiaoj@sina.com；通信作者：葛永斌(1975—)，男，寧夏吳忠人，教授，博士生導師，主要從事偏微分方程數(shù)值解法研究，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.