楊曉佳，田芳

(寧夏大學(xué) 數(shù)學(xué)統(tǒng)計(jì)學(xué)院,寧夏 銀川 750021)

?

一維非定常對(duì)流擴(kuò)散反應(yīng)方程的高精度緊致差分格式

楊曉佳，田芳

(寧夏大學(xué) 數(shù)學(xué)統(tǒng)計(jì)學(xué)院,寧夏 銀川 750021)

針對(duì)一維非定常對(duì)流擴(kuò)散反應(yīng)方程，首先推導(dǎo)了一種新的2層高精度緊致差分隱格式，其截?cái)嗾`差為O(τ2+τh2+h4)，即當(dāng)τ=O(h2)時(shí)，格式空間具有四階精度；然后采用Fourier分析方法分析了格式的穩(wěn)定性；最后通過(guò)數(shù)值算例驗(yàn)證了本文格式的精確性和可靠性.

對(duì)流擴(kuò)散反應(yīng)方程；非定常；緊致差分格式；隱式格式；高精度

MSC 2010：O 241.82

對(duì)流擴(kuò)散反應(yīng)方程是描述黏性流體流動(dòng)的非線性方程的線性化模型方程，該方程可用于描述許多自然現(xiàn)象，例如水和大氣中污染物質(zhì)濃度的擴(kuò)散，海水鹽度﹑溫度的擴(kuò)散，流體流動(dòng)和流體傳熱等過(guò)程，因此尋求此類(lèi)方程穩(wěn)定、實(shí)用的數(shù)值方法有著重要的理論與實(shí)際意義.目前大多數(shù)的文獻(xiàn)報(bào)道都針對(duì)的是對(duì)流擴(kuò)散方程模型，如葛永斌[1]等人通過(guò)時(shí)間導(dǎo)數(shù)項(xiàng)采用二階的向后歐拉差分公式，空間二階導(dǎo)數(shù)采用四階緊致差分公式，最終得到了一種求解含源項(xiàng)非定常變系數(shù)對(duì)流擴(kuò)散方程的3層四階緊致隱格式；Ding和Zhang[2]針對(duì)一維對(duì)流擴(kuò)散方程，采用半離散的方法和Padé逼近法，推導(dǎo)了一種截?cái)嗾`差為O(τ5+h4)的無(wú)條件穩(wěn)定的 隱格式.Tian和Yu[3]采用四階緊致指數(shù)差分公式和(2,2)Padé公式分別對(duì)空間變量和時(shí)間變量進(jìn)行離散，推導(dǎo)了一種求解一維非定常對(duì)流擴(kuò)散方程的四階精度的無(wú)條件穩(wěn)定的緊致差分隱格式；趙飛[4]等人針對(duì)一維非定常對(duì)流擴(kuò)散方程，考慮方程在n+1/2時(shí)刻的值，時(shí)間上利用中心差分格式，空間上采用截?cái)嗾`差余項(xiàng)修正的方法及加權(quán)平均的思想，最終推導(dǎo)出了非均勻網(wǎng)格上的一種2層高精度全隱式緊致差分格式；Piao等[5]針對(duì)一維非定常對(duì)流擴(kuò)散方程，采用3點(diǎn)緊致差分格式和快速單對(duì)角隱式Runge-Kutta方法分別對(duì)空間和時(shí)間進(jìn)行離散，得到了截?cái)嗾`差為O(τ4+h4)的無(wú)條件穩(wěn)定的隱格式.也有一些科研工作者針對(duì)對(duì)流擴(kuò)散反應(yīng)方程模型展開(kāi)研究，如Wo等[6]針對(duì)一維非定常對(duì)流擴(kuò)散反應(yīng)方程，采用Godunov方法，推導(dǎo)出了一種高精度的顯式差分格式.魏劍英[7]提出了一種求解非定常對(duì)流擴(kuò)散反應(yīng)方程的無(wú)條件穩(wěn)定的隱式差分格式，此格式的時(shí)間和空間方向均具有二階精度.Liao[8]先將一維非定常對(duì)流擴(kuò)散方程轉(zhuǎn)化為擴(kuò)散方程，然后采用四階緊致差分格式，最終推導(dǎo)出了一種截?cái)嗾`差為O(τ4+h4)的無(wú)條件穩(wěn)定的隱格式，該論文利用此格式計(jì)算了對(duì)流擴(kuò)散反應(yīng)問(wèn)題.趙秉新[9]針對(duì)一維對(duì)流擴(kuò)散反應(yīng)方程，通過(guò)指數(shù)變換將原方程變換為對(duì)流擴(kuò)散方程，對(duì)變換后的方程中的對(duì)流項(xiàng)和擴(kuò)散項(xiàng)分別采用高階迎風(fēng)緊致差分格式和對(duì)稱緊致格式進(jìn)行離散，在時(shí)間上采用四階Runge-Kutta方法進(jìn)行推進(jìn)，從而得到一種截?cái)嗾`差為O(τ4+h4)的緊致差分格式.楊錄峰等[10]針對(duì)對(duì)流擴(kuò)散反應(yīng)方程，采用消除對(duì)流項(xiàng)的處理技巧，結(jié)合四階Padé公式推導(dǎo)出了一種無(wú)條件穩(wěn)定的高精度差分格式.

本文針對(duì)一維非定常變系數(shù)對(duì)流擴(kuò)散反應(yīng)方程，首先采用泰勒級(jí)數(shù)展開(kāi)法和余項(xiàng)修正法，得到一種高精度緊致隱格式.然后采用Fourier分析法分析了格式的穩(wěn)定性.最后通過(guò)數(shù)值實(shí)驗(yàn)對(duì)本文方法的精確性和可靠性進(jìn)行驗(yàn)證.

1 高階緊致差分格式的建立

考慮如下一維變系數(shù)非定常對(duì)流擴(kuò)散反應(yīng)方程

(1)

其中，a為正常數(shù)，p(x,t)為對(duì)流項(xiàng)系數(shù)，q(x,t)為反應(yīng)項(xiàng)系數(shù)，f(x,t)為源項(xiàng).以τ表示時(shí)間步長(zhǎng)，空間取等間距網(wǎng)格，步長(zhǎng)用h表示.網(wǎng)格點(diǎn)為(xi,tn)，xi=ih,tn=nτ，i=0，1，…,N,n≥0.

考慮方程(1)在n第時(shí)刻的值并對(duì)時(shí)間和空間導(dǎo)數(shù)項(xiàng)分別采用向前和中心差分離散可得

(2)

).

(3)

為了得到時(shí)間二階和空間四階精度的緊致差分格式，需對(duì)式(3)中的時(shí)間二階導(dǎo)數(shù)項(xiàng)、空間三階和四階導(dǎo)數(shù)項(xiàng)進(jìn)行處理，為此對(duì)原方程(1)兩邊關(guān)于和分別求導(dǎo)并整理得

(4)

(5)

(6)

在式(4)、(5)和(6)中，全部變量關(guān)于空間的一階和二階導(dǎo)數(shù)項(xiàng)采用中心差分公式離散，未知量關(guān)于時(shí)間的一階導(dǎo)數(shù)項(xiàng)采用向前差分公式離散，各項(xiàng)系數(shù)p,q及右端項(xiàng)f關(guān)于時(shí)間的一階導(dǎo)數(shù)項(xiàng)采用二階向后歐拉差分公式離散，并整理得

(7)

(8)

(9)

(10)

將式(7)、(8)和(9)代入式(3)中并整理，可得到方程(1)的高階緊致差分格式

(11)

其中，

(12)

由推導(dǎo)過(guò)程可知，差分格式(11)的截?cái)嗾`差為O(τ2+τh2+h4)，即當(dāng)τ=h2時(shí)，格式在空間具有四階精度.注意到差分格式(11)對(duì)未知量u只涉及到了2層，但是對(duì)各項(xiàng)系數(shù)p,q和右端源項(xiàng)f涉及到了3層，由于這些函數(shù)都是已知的，因此可直接計(jì)算，不需要另外構(gòu)造格式啟動(dòng).

2 穩(wěn)定性分析

下面對(duì)格式(11)的穩(wěn)定性進(jìn)行分析.

考慮齊次的非定常對(duì)流擴(kuò)散反應(yīng)方程，即令方程(1)式右端項(xiàng)f=0，變形整理可得

(13)

取p(x,t)/a的上確界為Q，q(x,t)/a的上確界為C，令K=1/a，則(13)式變?yōu)?/p>

(14)

根據(jù)格式式(11)的推導(dǎo)過(guò)程，可得式(14)的離散格式為

(15)

其中，

(16)

(17)

則

(18)

其中，

A=(wi-1+wi+1)cosr+wi,B=(wi-1-wi+1)sinr,

F1=(gi-1+gi+1)cosr+gi,F2=(gi-1-gi+1)sinr.

在r-v平面內(nèi)作出|G|的等值線，圖1-4中分別給出了Pe=1、10、100、1 000，均取R=1時(shí)，|G|的等值線，|G|<1的區(qū)域表示格式穩(wěn)定，|G|>1的區(qū)域表示格式不穩(wěn)定.從圖中可以看出，差分格式式(11)對(duì)小波數(shù)以及大部分中小波數(shù)是穩(wěn)定的，圖3和圖4表明，隨著Pe數(shù)的增加，差分格式式(11)的穩(wěn)定區(qū)域不再發(fā)生變化.

3 數(shù)值算例

為了驗(yàn)證本文方法的精確性和可靠性，下面選取以下4個(gè)算例進(jìn)行數(shù)值實(shí)驗(yàn)，算例中的右端項(xiàng)f(x,t)及初邊值條件均由精確解給出， 并與已有文獻(xiàn)計(jì)算結(jié)果進(jìn)行比較.計(jì)算是用Fortran 77語(yǔ)言進(jìn)行編程，且在PC機(jī)上雙精度制下進(jìn)行的，離散后在每一個(gè)時(shí)間層上所得到的代數(shù)方程組可采用追趕法進(jìn)行求解.

圖1 當(dāng)R=1,Pe=1時(shí)，在r-v平面內(nèi)|G| 的等值線

圖2 當(dāng)R=1,Pe=10時(shí)，在r-v平面內(nèi)|G| 的等值線

圖3 當(dāng)R=1,Pe=100時(shí)，在r-v平面內(nèi)|G| 的等值線

圖4 當(dāng)R=1,Pe=1 000時(shí)，在r-v平面內(nèi)|G| 的等值線 Fig.4 Contours of the amplification factor |G| in the r-v plane for R=1,Pe=100

其中Ui表示點(diǎn)xi處的精確解，ui表示點(diǎn)xi處的數(shù)值解，Error(N1)和Error(N2)分別表示網(wǎng)格數(shù)為N1和N2時(shí)對(duì)應(yīng)的最大絕對(duì)誤差.

算例1[1]

其精確解為u(x,t)=e-tsin(x).

表1 問(wèn)題1當(dāng)τ=h2,t=0.5時(shí)的最大絕對(duì)誤差及收斂階

算例2[2,9]

其精確解為u(x,t)=e-tsin(x).

表2 問(wèn)題2當(dāng)τ=0.001,t=20時(shí)的L2誤差及收斂階

算例3

其精確解為u(x,t)=e5x+t(0.25+0.1π2)sinπx.

表3 問(wèn)題3當(dāng)τ=h2,t=0.5時(shí)的最大絕對(duì)誤差及收斂階

表4 問(wèn)題3當(dāng)N=32,t=1時(shí)對(duì)不同的網(wǎng)格比的最大絕對(duì)誤差

算例4

表5 問(wèn)題4當(dāng)τ=h2,t=1時(shí)的最大絕對(duì)誤差及收斂階

算例1是對(duì)流擴(kuò)散方程，表1給出了算例1當(dāng)τ=h2,t=0.5時(shí)的最大絕對(duì)誤差和收斂階的比較.從表中的計(jì)算結(jié)果可以看出，古典隱格式和C-N格式空間只具有二階精度，本文格式(11)空間具有四階精度，這與格式的理論分析相吻合.文獻(xiàn)[1]對(duì)此問(wèn)題空間也具有四階精度，但本文格式的計(jì)算結(jié)果比文獻(xiàn)[1]的計(jì)算結(jié)果更加精確，并且文獻(xiàn)[1]中的格式為3層格式，需要另外構(gòu)造格式來(lái)啟動(dòng)，而本文格式為2層格式，在計(jì)算過(guò)程中不需要另外構(gòu)造格式來(lái)啟動(dòng)，這使得本文格式使用時(shí)更加簡(jiǎn)便.算例2是對(duì)流擴(kuò)散方程，表2給出了算例2當(dāng)τ=0.001,t=20時(shí)的L2誤差和收斂階的比較.從表中的計(jì)算結(jié)果可以看出，本文格式與文獻(xiàn)[1,2,9]中的格式空間均具有四階精度，但文獻(xiàn)[1]中的格式是一個(gè)3層格式，需要另外構(gòu)造格式來(lái)啟動(dòng)，文獻(xiàn)[2]中的格式在求解過(guò)程中需要進(jìn)行矩陣求逆和乘積等運(yùn)算，當(dāng)網(wǎng)格數(shù)很大時(shí)，這2種格式的求解效率會(huì)明顯下降，本文格式的L2誤差比文獻(xiàn)[1,2,9]中格式的L2誤差小，這說(shuō)明本文格式具有更高的計(jì)算精度.算例3是一個(gè)對(duì)流擴(kuò)散反應(yīng)方程，由于文獻(xiàn)[1,2]中未給出關(guān)于此類(lèi)方程的計(jì)算格式，所以在此處未與文獻(xiàn)[1,2]做比較而僅與古典隱格式和C-N格式進(jìn)行了比較.表3和表4分別給出了算例3的計(jì)算結(jié)果，其中表3給出了當(dāng)τ=h2,t=0.5時(shí)，算例3的最大絕對(duì)誤差和收斂階的比較.從表中的計(jì)算結(jié)果可以看出，本文格式(11)對(duì)此問(wèn)題空間仍具有四階精度，而古典隱格式和C-N格式空間只具有二階精度，并且本文格式的最大絕對(duì)誤差比古典隱格式小1～4個(gè)數(shù)量級(jí)，比C-N格式小1～3個(gè)數(shù)量級(jí).表4給出了當(dāng)t=1時(shí)，網(wǎng)格數(shù)N取32，網(wǎng)格比λ分別取0.4,0.8,1.6,3.2,6.4時(shí)，算例3的最大絕對(duì)誤差，表中的計(jì)算結(jié)果與穩(wěn)定性的理論分析是一致的.算例4也是一個(gè)對(duì)流擴(kuò)散反應(yīng)問(wèn)題，表5給出了當(dāng)τ=h2,t=1時(shí)，算例4的最大絕對(duì)誤差和收斂階的比較.從表中的計(jì)算結(jié)果可以看出，本文格式對(duì)此問(wèn)題空間具有四階精度，而古典隱格式和C-N格式對(duì)此問(wèn)題的空間只具有二階精度，并且本文格式的最大絕對(duì)誤差比古典隱格式小2～4個(gè)數(shù)量級(jí)，比C-N格式小1～4個(gè)數(shù)量級(jí)，當(dāng)網(wǎng)格數(shù)較少時(shí)，古典隱格式和C-N格式的最大絕對(duì)誤差較大.

4 結(jié)論

首先建立了一種數(shù)值求解一維非定常對(duì)流擴(kuò)散反應(yīng)方程的2層全隱式緊致差分格式，格式在每個(gè)時(shí)間層上只用到了3個(gè)網(wǎng)格點(diǎn)，所以差分離散得到的代數(shù)方程組是三對(duì)角方程組，故可采用追趕法進(jìn)行求解.然后，通過(guò)Fourier分析方法分析了格式的穩(wěn)定性.最后選取了4個(gè)典型數(shù)值算例進(jìn)行了數(shù)值實(shí)驗(yàn).數(shù)值實(shí)驗(yàn)結(jié)果表明，當(dāng)τ=O(h2)時(shí)，本文格式在空間上可以達(dá)到4階精度，與文獻(xiàn)中的計(jì)算結(jié)果的比較顯示，本文方法總體上具有更高的計(jì)算精度.

此外，本文雖然只給出了一維對(duì)流擴(kuò)散反應(yīng)方程的高精度緊致隱格式，但此方法可以推廣到高維問(wèn)題的求解.對(duì)于高維方程的隱格式的求解，所得到的方程組不再是三對(duì)角線型的，所以不能直接采用追趕法，一般需要迭代求解.然而傳統(tǒng)的迭代法收斂速度很慢，為此可采用多重網(wǎng)格方法來(lái)加速收斂.文獻(xiàn)[11,12]給出了高維對(duì)流擴(kuò)散方程的高精度隱式差分格式及其多重網(wǎng)格算法，取得了較好的計(jì)算效果.從理論上講，將本文方法推廣到高維的對(duì)流擴(kuò)散反應(yīng)方程后，同樣可以采用多重網(wǎng)格方法進(jìn)行加速求解.

[1] 葛永斌,田振夫,吳文權(quán).含源項(xiàng)非定常對(duì)流擴(kuò)散方程的高精度緊致隱式差分方法[J].A輯,水動(dòng)力研究與進(jìn)展,2006,21(5):619-625. GE Y B,TIAN Z F,WU W Q.A high-order compact implic it difference method for the unsteady convection-diffusion equation with source term[J].Journal of Hydrodynamics,2006,21 (5):619-625.DOI:10.3969/j.issn.1000-4874.2006.05.010.

[2] DING H,ZHANG Y.A new difference scheme with high accuracy and absolute stability for solving convection-diffusion equatons[J].Comput Appl Math,2009,230:600-606.DOI.org/10.1016/j.cam.2008.12.015.

[3] TIAN Z F,YU P X.A high-order exponential scheme for solving 1D unsteady convection-diffusion Equations[J].J Comput Appl Math,2011,235:2477-2491.DOI.org/10.1016/j.cam.2010.11.001.

[4] 趙飛,陳建華,葛永斌.一維非定常對(duì)流擴(kuò)散方程非均勻網(wǎng)格上的高階緊致差分格式[J].西安理工大學(xué)學(xué)報(bào),2013,29 (4):475-480. ZHAO F,CHEN J H,GE Y B.A High order compact difference scheme on non-uniform grids for the 1D unsteady convection diffusion equation[J].Journal of Xi`an University of Technology,2013,29 (4):475-480.

[5] PIAO X,CHOI H J,KIM S D,et al.A fast singly diagonally implicit Runge-Kutta method for solving 1D unsteady convection-diffusion equations[J].Numer Methods Partial Differential Eq,2013,30:788-812.DOI.org/10.1002/num.21832.

[6] WO S,CHEN B M,WANG J.A high-order Godunov method for one-dimensional convection-diffusion-reaction problems[J].Numer Methods Partial Differential Eq,2000,16:495-512.DOI.org/10.1002/1098-2426(200011)16:6<495:aid-num1>3.0.co;2-s.

[7] 魏劍英.求解對(duì)流擴(kuò)散反應(yīng)方程的一種隱式差分格式[J].四川理工學(xué)院學(xué)報(bào)(自然科學(xué)版),2011,24(5):580-582. WEI J Y.An implicit scheme of the 1D convection-diffusion-reaction equation[J].Journal of Sichuan University of Science & Engineering(Natural Science Edition),2011,24(5):580-582.

[8] LIAO W.A compact high-order finite difference method for unsteady convection-diffusion equation[J].Intern J Comput Meth Eng Sci Mech,2012,13:135-145.DOI.org/10.1080/15502287.2012.660227.

[9] 趙秉新.一種求解一維對(duì)流擴(kuò)散反應(yīng)方程的高階緊致格式[J].重慶理工大學(xué)學(xué)報(bào)(自然科學(xué)版),2012,26(7),100-104. ZHAO B X.A High-order compact difference scheme for solving 1D convection-diffusion-reaction equation[J].Journal of Chongqin University of Technology(Natural Science),2012,26(7),100-104.

[10] 楊錄峰,李春光.一種求解對(duì)流擴(kuò)散反應(yīng)方程的高精度緊致差分格式[J].寧夏大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,34(2):101-109. YANG L F,LI C G.A High-order compact finite difference scheme for solving the convection diffusion reaction equations[J].Journal of Ningxia University(Natural Science Edition),2013,34(2):101-109

[11] 葛永斌,吳文權(quán),田振夫.二維對(duì)流擴(kuò)散方程的高精度全隱式多重網(wǎng)格方法[ J ].計(jì)算力學(xué)學(xué)報(bào),2004,21 (6):678-682. GE Y B,WU W Q,TIAN Z F.Multigrid method based on the high accuracy full implicit scheme of the convection diffusion equation[J].Chinese Journal of Computational Mechanics,2004,21 (6):678-682.DOI:10.3969/j.issn.1007-4708.2004.06.007.

[12] 葛永斌,田振夫,吳文權(quán).三維對(duì)流擴(kuò)散方程的高精度全隱式多重網(wǎng)格方法[ J ].計(jì)算力學(xué)學(xué)報(bào),2007,24 (2):181-186. GE Y B,TIAN Z F,WU W Q.Multigrid method based on high accuracy full implicit scheme of 3D convection diffusion equation[J].Chinese Journal of Computational Mechanics,2007,24 (2):181-186.DOI:10.3969/j.issn.1007-4708.2007.02.010.

(責(zé)任編輯：梁俊紅)

High order compact difference scheme for the one-dimensional unsteady convection diffusion reaction equation

YANG Xiaojia,TIAN Fang

(School of Mathematics and Statistis Science,Ningxia University,Yinchuan 750021,China)

A two-level high order compact finite difference implicit scheme is proposed to solve the one-dimensional unsteady convection diffusion reaction equation.The local truncation error of the scheme is O(τ2+τh2+h4),i.e.the scheme is the fourth order accuracy for space when τ=O(h2).Then,Fourier analysis method is used to prove the stability of the scheme.Finally,numerical experiments are conducted to verify the accuracy and the reliability of the present scheme.

convection diffusion reaction equation;unsteady;compact difference scheme;implicit scheme;high accuracy

10.3969/j.issn.1000-1565.2017.01.002

2016-03-28

國(guó)家自然科學(xué)基金資助項(xiàng)目(11361045;11161036);寧夏大學(xué)自然科學(xué)基金項(xiàng)目資助(ZR15014);寧夏大學(xué)研究生創(chuàng)新項(xiàng)目(GIP201620)

楊曉佳(1988—)，女，寧夏吳忠人，寧夏大學(xué)在讀碩士研究生，E-mail：yang_xiaoj@sina.com

田芳(1979—)，女，寧夏中寧人，寧夏大學(xué)副教授，主要從事偏微分方程數(shù)值解法的研究. E-mail：tianfsunny@126.com

O241.82

A

1000-1565(2017)01-0005-08