祁應(yīng)楠

(寧夏師范學(xué)院數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院，寧夏固原 756000)

?

求解一維拋物型方程的高精度有限差分方法

祁應(yīng)楠

(寧夏師范學(xué)院數(shù)學(xué)與計(jì)算機(jī)科學(xué)學(xué)院，寧夏固原 756000)

針對一維拋物型方程，采用樣條函數(shù)近似和Padé公式，構(gòu)造了一種高精度有限差分格式.該格式關(guān)于時(shí)間和空間均具有六階精度，并且從理論上被證明是無條件穩(wěn)定的.通過數(shù)值實(shí)驗(yàn)驗(yàn)證了本文方法的精確性和穩(wěn)定性，與文獻(xiàn)計(jì)算結(jié)果比較顯示,本文格式的計(jì)算結(jié)果更加精確.

拋物型方程；樣條函數(shù)；Padé逼近；高精度；有限差分法

拋物型方程是一類非常重要的偏微分方程，在擴(kuò)散、熱傳導(dǎo)、滲流等領(lǐng)域中有著廣泛的應(yīng)用，如絕熱過程中氣體通過多孔介質(zhì)的流動(dòng)問題、一個(gè)內(nèi)部有熱源的熱傳導(dǎo)過程、污染物濃度的擴(kuò)散問題等都可以用拋物型方程來描述.由于實(shí)際問題往往很難得到精確解，因此研究其數(shù)值解法具有重要的理論意義和實(shí)際應(yīng)用價(jià)值.

對拋物型方程很多情況下都采用有限差分法進(jìn)行求解.由于傳統(tǒng)的差分格式往往精度較低，一般僅具有一階或二階精度，因此不少研究者致力于發(fā)展高精度差分格式.馬明書等[1]針對一維拋物型方程建立了一族隱式差分格式，其時(shí)間為三階精度，空間為六階精度，格式是無條件穩(wěn)定的；Sun 和Zhang[2]針對一維拋物型方程，對空間導(dǎo)數(shù)項(xiàng)利用四階緊致差分，對時(shí)間項(xiàng)利用邊界值方法，構(gòu)造了無條件穩(wěn)定的四階精度差分方法.Sallam等[3]對空間變量利用中心差分公式進(jìn)行離散，得到關(guān)于時(shí)間變量的常微分方程組，然后基于C1三次樣條配置法得到了時(shí)間為四階精度而空間為二階精度的無條件穩(wěn)定的差分格式；葛永斌等[4]針對一維拋物型方程，利用Crank-Nicolson格式的思想和二階導(dǎo)數(shù)的四階緊致差分公式，分別建立了時(shí)間為二階、空間為四階和時(shí)間空間均為四階精度的緊致隱式差分格式，并利用離散Fourier分析法證明了前一種差分格式是無條件穩(wěn)定的，但后一種差分格式是無條件不穩(wěn)定的.Wen和Shao[5]采用區(qū)域分解法，構(gòu)造了一種求解一維和二維拋物型方程的高精度顯格式，其截?cái)嗾`差為O(τ3+h3)，并且是無條件穩(wěn)定的.Gao和Sun[6]建立了Neumann邊界條件下一維拋物型方程的一種緊致差分格式，此格式在內(nèi)網(wǎng)格點(diǎn)上的計(jì)算精度為時(shí)間二階和空間四階，而在邊界點(diǎn)上的計(jì)算精度為時(shí)間二階和空間三階.劉蕤和高銳敏[7]基于廣義梯形公式和四次樣條函數(shù)，建立了Neumann邊界條件下的一維拋物型方程的一族含參數(shù)的隱式差分格式，該格式在空間方向上達(dá)到了四階精度，但時(shí)間方向上僅有二階精度.當(dāng)取特殊參數(shù)時(shí)，該格式的時(shí)間精度可提高到三階.陳國玉等[8]建立了一維拋物型方程的一種條件穩(wěn)定的三層高精度差分格式，在時(shí)間方向上具有二階精度，而空間方向上具有四階精度.

本文針對一維拋物型方程，采用樣條函數(shù)近似和Padé公式，構(gòu)造一種時(shí)間和空間均具有六階精度的有限差分格式，并且從理論和數(shù)值實(shí)驗(yàn)兩方面分析和驗(yàn)證本文所提格式的精確性和穩(wěn)定性.考慮一維拋物型方程

(1)

初始條件

(2)

邊界條件

(3)

其中φ(x)是已知函數(shù)，a>0為擴(kuò)散項(xiàng)系數(shù)，α,β是常數(shù).

1 高階有限差分格式

(4)

(5)

(6)

(7)

根據(jù)文獻(xiàn)[9]的結(jié)論，有

(8)

(9)

組合(1),(8)和(9)式，即可得到

(10)

即

(11)

又因?yàn)?/p>

(12)

于是有

(13)

即

(14)

舍掉高階截?cái)嗾`差項(xiàng)，則方程(1)的半離散格式可寫為

(15)

(15)式包含N-3個(gè)方程，但根據(jù)計(jì)算域離散有N-1個(gè)未知節(jié)點(diǎn)，因此還需補(bǔ)充兩個(gè)邊界點(diǎn)方程，方程組才有唯一解.為了得到與(15)式相匹配的六階精度格式，定義以下恒等式[10]：

(16)

(17)

利用泰勒級數(shù)展開，可以得到以下邊界方程：

4u0-7u1+2u2+u3=

(18)

(19)

利用(18)和(19)式，可以得到問題(1)～(3)的兩個(gè)邊界方程分別為：

(20)

(21)

以上兩種格式的截?cái)嗾`差均為O(h6).

由(3)式得：u0(t)=α,uN(t)=β,(ut)0(t)=0,(ut)N(t)=0.令

結(jié)合(15),(20)和(21)式，有

(22)

其中

由于矩陣A是嚴(yán)格對角占優(yōu)的，B是不可約對角占優(yōu)矩陣，那么矩陣A和B都是非奇異的[11]，結(jié)合初始條件(2)，有

(23)

其中

則(23)式的解為

(24)

其中

I是(N-1)×(N-1)階單位矩陣,則(24)式可以變形為

(25)

令X(t)=U(t)-B-1b,X0=U0-B-1b,C=A-1B,則(25)式可重新寫為

(26)

令τ為時(shí)間步長，則(26)式在點(diǎn)(xi,tn)處的離散結(jié)果為

(27)

其中tn=nτ,n=0,1,2,….令X(tn+1)=Xn+1,X(tn)=Xn，則

(28)

其中n=0,1,2,….由于e-τ C是關(guān)于矩陣τC的一個(gè)無窮級數(shù)，為了逼近e-τ C，我們采用如下(3,3)Padé逼近[12]：

(29)

即采用

(120I-60τC+12τ2C2-τ3C3)×

(30)

來逼近e-τ C，那么(28)式的離散格式為

(31)

其中n=0,1,2,….從格式構(gòu)造過程易知格式(31)的截?cái)嗾`差為O(τ6+h6).

2 格式的穩(wěn)定性分析

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

D(120I-60τC+12τ2C2-τ3C3)×

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

記

Ej=120-60τ(aj+bji)+12τ2(aj+bji)2-

Fj=120+60τ(aj+bji)+12τ2(aj+bji)2+

則有

于是有

由于aj>0，因此有

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

3 數(shù)值算例

為了驗(yàn)證本文格式(31)的性能，我們選用如下幾個(gè)問題進(jìn)行數(shù)值實(shí)驗(yàn).收斂階、最大和平均絕對誤差的定義如下：

其中，Max_erri(i=1,2)是空間網(wǎng)格步長為hi(i=1,2)時(shí)的最大絕對誤差，Ave_erri(i=1,2)是空間網(wǎng)格步長是hi(i=1,2)時(shí)的平均絕對誤差，u(xj)是在xj點(diǎn)處的精確解，Uj是在xj點(diǎn)處的數(shù)值解.

問題1[2]ut=(1/π2)uxx,00,u(x,0)=sin(πx),u(0,t)=u(1,t)=0.該問題的精確解為u(x,t)=e-tsin(πx).

表1 問題1在τ=h,t=1.0時(shí)的最大絕對誤差及收斂階

表2 問題1在h=1/40,t=1.0時(shí)對不同網(wǎng)格比λ的最大絕對誤差

表1給出了空間步長h分別取1/8,1/16,1/32,1/64時(shí)，問題1在時(shí)間t=1.0時(shí)刻C-N格式、文獻(xiàn)[2]格式和本文格式的最大絕對誤差及收斂階.可以看出，在相同網(wǎng)格等分?jǐn)?shù)下，由本文格式得到的最大絕對誤差比C-N格式和文獻(xiàn)[2]格式小好幾個(gè)數(shù)量級.當(dāng)τ=O(h)時(shí)，本文格式達(dá)到了六階精度，而C-N格式為二階精度，文獻(xiàn)[2]格式為四階精度,且只有在τ=O(h2)時(shí)精度才能達(dá)到四階.

表2給出了問題1在不同網(wǎng)格比λ下的最大絕對誤差，以此來驗(yàn)證三種格式的穩(wěn)定性.可以看出，計(jì)算均是穩(wěn)定的,并且本文格式的計(jì)算結(jié)果比C-N格式和文獻(xiàn)[2]格式精確得多.另外，從理論上講本文方法時(shí)間和空間均具有六階精度，由于本文空間導(dǎo)數(shù)采用了樣條方法(本質(zhì)上是多項(xiàng)式逼近)，而時(shí)間導(dǎo)數(shù)采用了(3,3)Padé近似(本質(zhì)上是指數(shù)型逼近)，因此關(guān)于時(shí)間和空間的截?cái)嗾`差主項(xiàng)系數(shù)具有不同的函數(shù)特征，從而不能保證計(jì)算精度正好是六階，所以會出現(xiàn)一定的波動(dòng)現(xiàn)象.

問題2[3]ut=uxx,0

表3 問題2在λ=1/2,t=2.0時(shí)的平均絕對誤差及收斂階

表4 問題2在λ=1,t=2.0時(shí)的平均絕對誤差及收斂階

表3和表4分別給出了當(dāng)λ=τ/h2=1/2和1，空間步長h分別取π/8,π/16,π/32,π/64時(shí)，問題2在時(shí)間t=2.0時(shí)刻C-N格式、文獻(xiàn)[3]格式和本文格式的平均絕對誤差及收斂階.從表中可以看出，本文格式的結(jié)果比C-N格式和文獻(xiàn)[3]格式精確得多.本文格式達(dá)到了六階精度，而C-N格式和文獻(xiàn)[3]格式僅為二階精度.

問題3[4]ut=auxx,00,u(x,0)=sin(πx),u(0,t)=u(1,t)=0.該問題的精確解為u(x,t)=e-aπ2tsin(πx).

表5給出了問題3當(dāng)a=1,τ=h2,h分別取1/8,1/16,1/32,1/64時(shí)，在時(shí)間t=0.25時(shí)刻C-N格式、文獻(xiàn)[4]格式和本文格式的最大絕對誤差及收斂階.可以看出，C-N格式僅有二階精度，文獻(xiàn)[4]格式具有四階精度，而本文格式達(dá)到了六階精度.表6給出了問題3當(dāng)a=1,h=1/32時(shí)不同時(shí)間步長τ在t=0.5時(shí)刻的最大絕對誤差.可以看出，對該問題，網(wǎng)格比λ從0.8變化到51.2，三種格式的計(jì)算均是穩(wěn)定的,并且本文格式的計(jì)算結(jié)果比C-N格式和文獻(xiàn)[4]格式精確得多.

表5 問題3在a=1,τ=h2,t=0.25時(shí)的最大絕對誤差及收斂階

表6 問題3在a=1,h=1/32,t=0.5對不同網(wǎng)格比λ的最大絕對誤差

4 結(jié)論

本文針對一維拋物型方程，首先對空間二階導(dǎo)數(shù)采用六次樣條函數(shù)進(jìn)行逼近，從而得到了半離散格式；接下來利用常微分方程精確解公式，得到關(guān)于時(shí)間t的指數(shù)矩陣表達(dá)式；然后采用(3,3)Padé逼近得到了一種求解一維拋物型方程的六階精度有限差分格式，并從理論上證明了該格式是絕對穩(wěn)定的.最后，數(shù)值實(shí)驗(yàn)驗(yàn)證了本文格式在時(shí)間和空間上的精度都達(dá)到了六階，并且比傳統(tǒng)的C-N格式和文獻(xiàn)[3-5]格式的計(jì)算結(jié)果更加精確.

[1] 馬明書,李改弟.拋物型方程的一個(gè)新的高精度恒穩(wěn)定的隱式差分格式[J].數(shù)學(xué)的實(shí)踐與認(rèn)識,2001,31(3):355.

[2] UN H W,ZHANG J.A high-order compact boundary value method for solving one-dimensional heat equations[J].NumericalMethodsforPartialDifferentialEquations,2003,19:846.

[3] SALLAM S,NAIM A M,ABDEL-AZIZ M R.Unconditionally stableC1-cubic spline collocation method for solving parabolic equations[J].InternationalJournalofComputerMathematics,2004,81(7):813.

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

[5] WEN R H,SHAO H Z.Domain decomposition schemes with high-order accuracy and unconditional stability[J].AppliedMathematicsandComputation,2013,219:6170.

[6] GAO G H,SUN Z Z.Compact difference schemes for heat equation with Neumann boundary conditions(Ⅱ)[J].NumericalMethodsforPartialDifferentialEquations,2013,29:1459.

[7] 劉蕤,高銳敏.帶Neumann邊界條件的拋物型方程的樣條差分方法[J].鄭州大學(xué)學(xué)報(bào)(理學(xué)版),2013,45(3):37.

[8] 陳國玉,謝英超,程燕.熱傳導(dǎo)方程的一個(gè)三層差分格式[J].四川兵工學(xué)報(bào),2014,35(7):143.

[9] KHAN A,KHAN I,AZIZ T.Sextic spline solution of a singularly perturbed boundary-value problems[J].AppliedMathematicsandComputation,2006,181:432.

[10] RASHIDINIA J,JALILIAN R,MOHAMMADI R,et al.Sextic spline method for the solution of a system of obstacle problems[J].AppliedMathematicsandComputation,2007,190:1669.

[11] 孫文瑜,杜其奎,陳金如.計(jì)算方法[M].北京:科學(xué)出版社,2007.

[13] 金升平,熊方方,李瓊.矩陣的實(shí)特征值為正的條件與判斷[J].重慶理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2015,25(1):117.

(責(zé)任編輯 馬宇鴻)

A high order finite difference method for solving the one-dimensional parabolic equation

QI Ying-nan

(College of Mathematics and Computer Science,Ningxia Normal University,Guyuan 756000,Ningxia,China)

A high order finite difference scheme for solving the one-dimensional parabolic equation is developed based on the Sextic spline functions for the spatial derivative and (3,3) Padé approximation for the temporal derivative.The scheme is sixth order accuracy in both time and space,and the performance of unconditional stability is proved in theory.The accuracy and effectiveness of the present method are validated by three numerical experiments,the results show that the present method is more accurate than the methods in the literature.

parabolic equation;spline function;Padé approximation;high order accuracy;finite difference method

10.16783/j.cnki.nwnuz.2016.06.006

2015-12-03;修改稿收到日期:2016-06-03

寧夏高等學(xué)校科學(xué)研究資助項(xiàng)目(NGY2015115);寧夏自然科學(xué)基金資助項(xiàng)目(NZ15259，NZ16251);寧夏師范學(xué)院科研項(xiàng)目(NXSFZD1603)

祁應(yīng)楠(1978—)，男，寧夏固原人，副教授，碩士.主要研究方向?yàn)槠⒎址匠虜?shù)值解法、計(jì)算流體力學(xué).

E-mail:421270380@qq.com

O 241.82

A

1001-988Ⅹ(2016)06-0029-06