楊曉佳，　葛永斌

(寧夏大學(xué) 數(shù)學(xué)計(jì)算機(jī)學(xué)院　寧夏 銀川 750021)

?

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

楊曉佳，葛永斌

(寧夏大學(xué) 數(shù)學(xué)計(jì)算機(jī)學(xué)院寧夏 銀川 750021)

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

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

0引言

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

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

1差分格式的建立

考慮一維含源項(xiàng)的擴(kuò)散方程的初邊值問題:

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)為擴(kuò)散項(xiàng)系數(shù)；f(x,t)為非齊次項(xiàng)；a(t),b(t),φ(x)均為已知函數(shù).

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

(4)

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

(5)

(6)

令

于是有

(7)

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

(8)

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

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

(9)

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

(10)

(11)

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

(12)

(13)

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

(14)

式(14)即為所構(gòu)造的差分格式,通過對空間的離散過程不難發(fā)現(xiàn)，其具有四階精度. 另外，由文獻(xiàn)[1]的研究結(jié)論可知，時(shí)間具有二階精度，因此本文所構(gòu)造的差分格式(14)的截?cái)嗾`差為O(τ2+h4).

2穩(wěn)定性分析

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

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

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

所以xTSx>0， 而

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

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

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

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

記

則有

那么

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

3數(shù)值算例

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

Error1和Error2分別為空間網(wǎng)格步長為h1和h2時(shí)的最大絕對誤差.

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

表1　問題1，當(dāng)τ=h2,t=1時(shí)對不同N的最大絕對誤差及收斂階

表2　問題1，當(dāng)N=32，t=0.5時(shí)對不同的r最大絕對誤差

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

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

表5　問題4，當(dāng)N=32,τ=h2，t=5時(shí)的數(shù)值解和精確解

表6　問題4 當(dāng)N=32，t=1時(shí)對不同r的最大絕對誤差

4結(jié)論

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

參考文獻(xià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]. 計(jì)算數(shù)學(xué),1982,4(2): 204—213.

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

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

[6]詹涌強(qiáng). 求解拋物型方程的一族六點(diǎn)隱式差分格式[J]. 安徽大學(xué)學(xué)報(bào)(自然科學(xué)版),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]開依沙爾熱合曼, 努爾買買提黑力力. 求解擴(kuò)散方程的二級四階隱式Runge-Kutta方法[J]. 湖北大學(xué)學(xué)報(bào)(自然科學(xué)版), 2014, 36(5): 476—480.

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

[10]周敏, 高學(xué)軍, 董超. 解拋物型方程的八點(diǎn)隱式差分格式[J]. 廣東工業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版), 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.

(責(zé)任編輯：方惠敏)

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

基金項(xiàng)目:國家自然科學(xué)基金資助項(xiàng)目(11361045, 11161036)；寧夏高等學(xué)?？茖W(xué)技術(shù)研究項(xiàng)目(NGY2013019).

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

中圖分類號：O241.82

文獻(xiàn)標(biāo)志碼：A

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

DOI:10.3969/j.issn.1671-6841.201506013

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