魏劍英，王燕，葛永斌

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

非穩(wěn)態(tài)對(duì)流擴(kuò)散方程是一類基本的發(fā)展方程，在生態(tài)環(huán)境、流體力學(xué)、生物數(shù)學(xué)等領(lǐng)域都有廣泛應(yīng)用．由于解析解通常很難求出，因此，研究各種有效實(shí)用的數(shù)值方法對(duì)此類方程進(jìn)行求解就顯得極為重要。

有關(guān)非穩(wěn)態(tài)對(duì)流擴(kuò)散方程的有限差分法的文獻(xiàn)報(bào)道已有很多，如文獻(xiàn)[1]利用半離散法和Padé逼近，推導(dǎo)了一維方程的截?cái)嗾`差為O(τ5+h4)的隱格式，且是無(wú)條件穩(wěn)定的； 文獻(xiàn)[2]利用Taylor級(jí)數(shù)展開(kāi)和待定系數(shù)法構(gòu)造了一致四階緊致格式并結(jié)合三階TVD-Runge-Kutta方法求解了一維方程； 文獻(xiàn)[3]針對(duì)二維方程，對(duì)空間導(dǎo)數(shù)項(xiàng)和時(shí)間導(dǎo)數(shù)項(xiàng)分別采用四階離散和加權(quán)平均的方法，推導(dǎo)出一種截?cái)嗾`差為O(τ2+h4)的格式； 文獻(xiàn)[4]采用半離散的方法，從一維定常問(wèn)題出發(fā)，利用四階緊致差分算子和Crank-Nicolson(C-N)格式，推導(dǎo)了一種求解二維常系數(shù)非穩(wěn)態(tài)對(duì)流擴(kuò)散方程的差分格式，截?cái)嗾`差為O(τ2+h4)，且是無(wú)條件穩(wěn)定的； 文獻(xiàn)[5]針對(duì)此類方程，構(gòu)造了有理型緊致交替方向隱式(ADI)格式，截?cái)嗾`差為O(τ2+h4)，且是無(wú)條件穩(wěn)定的； 文獻(xiàn)[6]還提出了該方程無(wú)條件穩(wěn)定的指數(shù)型緊致ADI差分格式，其截?cái)嗾`差依然為O(τ2+h4)； 文獻(xiàn)[7]利用指數(shù)變換消除對(duì)流項(xiàng)后轉(zhuǎn)化為擴(kuò)散方程，再利用緊致公式和擴(kuò)展的Simpson公式構(gòu)造了一維方程的高階格式，截?cái)嗾`差為O(τ4+h4)，且是無(wú)條件穩(wěn)定的； 文獻(xiàn)[8]針對(duì)二維非穩(wěn)態(tài)常系數(shù)方程，空間方向直接采用六階組合緊致差分公式進(jìn)行計(jì)算，時(shí)間方向用C-N格式離散，所提格式的截?cái)嗾`差為O(τ2+h6)，盡管該格式空間達(dá)到了六階精度，但是由于其時(shí)間只有二階精度，因此為了保證空間精度達(dá)到六階，其計(jì)算所需要采取的時(shí)間步長(zhǎng)必須為O(h3)，即必須采用較小的時(shí)間步長(zhǎng)，并且該方法僅適用于常系數(shù)問(wèn)題。

本文針對(duì)二維非穩(wěn)態(tài)變系數(shù)對(duì)流擴(kuò)散方程，首先對(duì)空間二階導(dǎo)數(shù)采用一階導(dǎo)數(shù)的四階逼近公式，而一階導(dǎo)數(shù)采用四階Padé逼近，時(shí)間項(xiàng)采用二階向后差分公式(BDF)，得到一個(gè)時(shí)間二階、空間四階精度的緊致差分格式； 然后利用Taylor級(jí)數(shù)展開(kāi)，對(duì)一階和二階空間導(dǎo)數(shù)項(xiàng)采用六階緊致差分公式，時(shí)間項(xiàng)采用三階BDF得到一個(gè)時(shí)間三階、空間六階精度的緊致差分格式．由于采用了向后差分，因此兩種格式均是無(wú)條件穩(wěn)定的。

1 差分格式的建立

考慮如下二維非穩(wěn)態(tài)對(duì)流擴(kuò)散方程：

ut-α(uxx+uyy)+p(x，y，t)ux+q(x，y，t)uy=f(x，y，t) (x，y)∈Ω，t≥0

(1)

初始條件：

u(x，y，0)=g(x，y) (x，y)∈Ω

(2)

邊界條件：

u(x，y，t)=s(x，y，t) (x，y)∈?Ω，t>0

(3)

其中：Ω∈[a1，a2]×[b1，b2]為R2上的矩形區(qū)域，?Ω為Ω的邊界，u(x，y，t)為待求未知量，p(x，y，t)，q(x，y，t)分別為x，y方向的對(duì)流項(xiàng)系數(shù)，α為擴(kuò)散項(xiàng)系數(shù)(常數(shù))，f(x，y，t)為源項(xiàng)，且p(x，y，t)，q(x，y，t)，f(x，y，t)，g(x，y)，s(x，y，t)均為已知函數(shù)。

(4)

(5)

1.1 四階精度緊致格式

將方程(1)改寫(xiě)為

ut=f(x，y，t)+α(uxx+uyy)-p(x，y，t)ux-q(x，y，t)uy

(6)

對(duì)二階導(dǎo)數(shù)uxx和uyy采用如下逼近：

(7)

(8)

將(7)，(8)式代入(6)式，考慮其在點(diǎn)(xi，yj)的值，得到

(9)

為了使(9)式空間保持四階精度，對(duì)ux和uy的計(jì)算采用如下四階Padé公式[9]

(10)

(11)

考慮(9)式在第n+1時(shí)間層的離散值，對(duì)ut采用如下二階BDF[10-11]進(jìn)行離散

(12)

舍去高階項(xiàng)，可得

(13)

(14)

其中G為(9)式的右端項(xiàng)。

在文獻(xiàn)[12]中已經(jīng)證明當(dāng)K≤6(K代表精度)時(shí)，BDF是無(wú)條件穩(wěn)定的．另外，由式(13)可以看出，在每個(gè)時(shí)間步上，該格式為5點(diǎn)模板的全隱格式．為了得到u(x，y，t)的計(jì)算值，(14)式中出現(xiàn)的ux和uy除需知道其內(nèi)點(diǎn)值，還需知道其在邊界點(diǎn)處的值，為保持與內(nèi)點(diǎn)差分格式同樣的精度，對(duì)邊界上的ux，uy采用一致四階邊界條件[2]計(jì)算。

1.2 六階精度緊致格式

為了得到方程(1)的六階精度的格式，需由Taylor級(jí)數(shù)展開(kāi)得到一、二階導(dǎo)數(shù)的六階近似公式

(15)

(16)

(17)

(18)

將(15)-(18)式代入(6)式，有

(19)

(20)

(21)

(22)

(23)

(24)

(25)

將(20)-(25)式代入(19)式，考慮其在點(diǎn)(xi，yj)的值，即可得到如下半離散的緊致格式

(26)

為了使(26)式具有六階精度，對(duì)ux和uy采用文獻(xiàn)[9]中的六階差分公式計(jì)算，uxx和uyy采用文獻(xiàn)[13]中的六階差分公式計(jì)算

(27)

(28)

(29)

(30)

考慮(26)式在第n+1時(shí)間層的值，對(duì)ut采用三階BDF[10]進(jìn)行離散

(31)

略去高階項(xiàng)，可得

(32)

2 數(shù)值實(shí)驗(yàn)

本文將利用以下4個(gè)問(wèn)題進(jìn)行數(shù)值實(shí)驗(yàn)，驗(yàn)證(2，4)格式和(3，6)格式的穩(wěn)定性和有效性，文中涉及到的誤差及收斂階定義如下：

1) 最大絕對(duì)誤差：

2)L2范數(shù)誤差：

3) 收斂階：

問(wèn)題1考慮如下非齊次對(duì)流擴(kuò)散問(wèn)題

其精確解為u(x，t)=e-txy(1-x)(1-y)，初邊值條件及右端項(xiàng)f(x，y，t)均由精確解給出。

表1給出了問(wèn)題1當(dāng)τ=h2，t=0.5時(shí)的計(jì)算結(jié)果．從表中可以看出，本文(2，4)格式和文獻(xiàn)[6]格式在空間上均達(dá)到了四階精度，但本文(2，4)格式的計(jì)算結(jié)果比BTCS格式和文獻(xiàn)[6]格式更好； 本文(3，6)格式在空間上達(dá)到了六階精度(網(wǎng)格為64時(shí)的數(shù)值已經(jīng)達(dá)到機(jī)器精度，所以收斂階受到影響不準(zhǔn)確)，且最大絕對(duì)誤差比(2，4)格式的小2～3個(gè)數(shù)量級(jí)，計(jì)算結(jié)果更精確．當(dāng)h=0.01，t=0.5時(shí)，不同時(shí)間步長(zhǎng)τ下的L∞誤差、收斂階及CPU時(shí)間由表2給出．從表中可以看出，本文(2，4)格式時(shí)間上達(dá)到的精度為二階，本文(3，6)格式時(shí)間上達(dá)到的精度為三階，這與理論分析一致； 另外當(dāng)τ較大時(shí)計(jì)算收斂慢，所用的CPU時(shí)間相對(duì)較長(zhǎng)，計(jì)算過(guò)程中(3，6)格式需要計(jì)算在節(jié)點(diǎn)處的一、二階導(dǎo)數(shù)值，因此所用的計(jì)算時(shí)間比(2，4)格式長(zhǎng)。

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

表2 問(wèn)題1當(dāng)h=0.01，t=0.5時(shí)的最大絕對(duì)誤差L∞，收斂階及CPU時(shí)間

問(wèn)題2考慮如下變系數(shù)對(duì)流擴(kuò)散問(wèn)題

其精確解為u(x，t)=txy(1-x)(1-y)ex+y，初邊值條件及右端項(xiàng)f(x，y，t)均由精確解給出。

表3給出了問(wèn)題2當(dāng)τ=h2，t=0.25時(shí)本文格式與C-N格式和BTCS格式的L∞誤差及收斂階的比較．從表中可以看出，C-N格式和BTCS格式在空間上達(dá)到的精度均為二階，本文(2，4)格式在空間上達(dá)到的精度為四階，而本文(3，6)格式空間上達(dá)到的精度為六階，且(3，6)格式的計(jì)算結(jié)果誤差更小，充分驗(yàn)證了本文兩種格式的精確性和有效性．另外，上述兩個(gè)例子的計(jì)算結(jié)果進(jìn)一步說(shuō)明低精度的啟動(dòng)步并不影響格式的整體精度。

表3 問(wèn)題2當(dāng)τ=h2，t=0.25時(shí)的L∞誤差及收斂階

問(wèn)題3考慮如下高斯脈沖問(wèn)題

圖1 t=1，τ=2.5×10-3，Pe=2，1.2≤x，y≤1.8時(shí)的等值線圖

圖2 t=0.001，τ=2.5×10-6，Pe=2000，1.2≤x，y≤1.8時(shí)的等值線圖

表4 問(wèn)題3當(dāng)h=0.02時(shí)在不同參數(shù)下的L∞誤差和L2誤差

問(wèn)題4考慮如下非線性方程

表5給出了問(wèn)題4當(dāng)t=1，τ=h2時(shí)，本文(2，4)格式和(3，6)格式在α分別取1和0.1時(shí)的最大絕對(duì)誤差與收斂階．從表中可以看出，本文(2，4)格式在空間上達(dá)到的精度是四階，(3，6)格式達(dá)到的精度是六階(網(wǎng)格為64，α=1時(shí)最大絕對(duì)誤差的值已經(jīng)達(dá)到機(jī)器精度)，這與理論推導(dǎo)是一致的，從表中還可以看出，α=1時(shí)所用的CPU時(shí)間比α=0.1時(shí)的長(zhǎng)，(3，6)格式的CPU時(shí)間比(2，4)格式的長(zhǎng)．表6給出了當(dāng)N=7，τ=0.01，α=1時(shí)，t=15和t=20的計(jì)算解與精確解的絕對(duì)誤差．表中數(shù)據(jù)表明，本文(2，4)格式的計(jì)算誤差比文獻(xiàn)[14]格式小5～6個(gè)數(shù)量級(jí)，與文獻(xiàn)[15]格式的計(jì)算結(jié)果大致相當(dāng)； 本文(3，6)格式的計(jì)算結(jié)果比文獻(xiàn)[14]格式小7～9個(gè)數(shù)量級(jí)，比文獻(xiàn)[15]格式小2～3個(gè)數(shù)量級(jí)．因此，對(duì)于此類非齊次邊界的Burgers方程，本文兩種高精度緊致差分格式非常有效．圖3給出了N=20，t=1，τ=0.01，α=0.1時(shí)本文兩種格式的絕對(duì)誤差和文獻(xiàn)[15]格式的絕對(duì)誤差比較，從圖中可以看出，本文(2，4)格式的絕對(duì)誤差與文獻(xiàn)[15]格式的絕對(duì)誤差大致相當(dāng)，但本文(3，6)格式的計(jì)算結(jié)果更加精確，絕對(duì)誤差更小。

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

表6 問(wèn)題4當(dāng)N=7，τ=0.01，α=1時(shí)的絕對(duì)誤差

續(xù)表6

圖3 N=20，t=1，τ=0.01，α=0.1時(shí)的精確解與絕對(duì)誤差

3 結(jié)論