侯波,葛永斌

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

1.引言

對流方程在生物數(shù)學(xué)、能源開發(fā)、空氣動力學(xué)等許多領(lǐng)域都具有十分廣泛的應(yīng)用,因此求解該類方程具有非常重要的理論價值和實際意義.然而,由于實際問題通常十分復(fù)雜,往往難以求得精確解,因此研究其精確穩(wěn)定的數(shù)值解法是十分必要的.

針對對流方程國內(nèi)外很多學(xué)者提出了很多的數(shù)值方法.如張?zhí)斓潞蛯O傳灼[1]針對一維對流方程采用待定系數(shù)法,得到了兩層四點格式和四階六點格式,并且是無條件穩(wěn)定的,該方法適用于在點數(shù)確定的前提下,得到精度高的差分格式; 于志玲和朱少紅[2]針對一維問題建立了中間層為兩個節(jié)點的三層顯格式,其截斷誤差為O(τ2+h2); 曾文平[3]針對一維對流方程推導(dǎo)出了一種兩層半顯式格式，其截斷誤差為O(τ2+h2),該格式是無條件穩(wěn)定的.姚朝輝等人[5]將二階的迎風(fēng)格式和中心差分格式進行加權(quán)得到了WSUC格式,該格式是無條件穩(wěn)定的;但該格式時間方向和空間方向僅有二階精度.湯寒松等人[6]通過立方插值擬質(zhì)點方法(CIP 方法),給出了一些保單調(diào)的CIP格式; Erdogan[9]針對一維的對流方程推導(dǎo)出了一種指數(shù)擬合的差分格式,其截斷誤差為O(τ2+h2); Bourchtein[10]構(gòu)造了對流方程的三層五點中心型蛙跳格式,該格式的截斷誤差為O(τ4+h4); 即該格式時間和空間均具有四階精度,但是該格式是三層的,空間方向需要五個點,并且是條件穩(wěn)定的; Kim[11]構(gòu)造了多層無耗散的迎風(fēng)蛙跳格式,即時間和空間分別具有二階、四階、六階精度,但格式為三層甚至是四層的,并且六階格式空間方向最多需要五個點,給靠近邊界的內(nèi)點的計算帶來困難.

綜上所述,文獻中已經(jīng)有的數(shù)值計算方法大多為低階精度的,而高精度方法涉及多個時間層,需要一個或多個時間啟動步,或者空間方向的網(wǎng)格節(jié)點多于三個,這都給計算造成困難或不便.為此本文將構(gòu)造一種緊致格式,這里緊致格式的定義為對時間導(dǎo)數(shù)項的離散采用不超過三個時間層,而對空間導(dǎo)數(shù)項的離散采用不超過三個網(wǎng)格點,時間和空間即可以達到高階精度（三階及三階以上）的格式.本文擬構(gòu)造的格式時間方向僅用到兩個時間層上的函數(shù)值,在每個時間層上僅涉及到三個空間網(wǎng)格點,格式時間和空間具有整體的四階精度.該格式的優(yōu)點是無須啟動步的計算,并且在對靠近邊界點的計算時,不會用到計算域以外的網(wǎng)格節(jié)點.此外該格式為無條件穩(wěn)定的,可以采用比較大的時間步長進行計算.最后通過數(shù)值實驗驗證本文格式的精確性和穩(wěn)定性.

2.差分格式的建立

考慮如下一維對流方程：

給定初始條件為:

給定周期性邊界條件為:

其中,u(x,t)為未知函數(shù),f為非齊次項,a為對流項系數(shù),φ(x)為已知函數(shù).

將求解區(qū)域[b,c]等距剖分為N個子區(qū)間：b=x0,x1,··· ,xN?1,xN=c,并且定義時間也采用等距剖分,步長用τ表示.在本文中,我們利用分別表示u在(xi,tn),點處的函數(shù)值.

其中

同理可得:

將(2.7)和(2.8)代入(2.4)整理可得:

為了使該格式在時間方向和空間方向上均達到四階精度,須對(2.9)式中的項進行二階的離散,同時為了保證本文格式的緊致性,即空間方向不超過三個網(wǎng)格點,我們對(2.1)式進行如下變形:

從而可得:

將(2.13)代入(2.11)得:

即,

即,

由推導(dǎo)過程可知,該格式的截斷誤差為O(τ4+τ2h2+h4),即格式(2.15) 在時間和空間上均可達到四階精度.我們注意到,格式為兩層格式,并且格式每層僅用到三個網(wǎng)格點,形成的代數(shù)方程組系數(shù)矩陣為循環(huán)三對角矩陣,可采用追趕法進行求解[8],同時由于要求未知時間層上(第n+1層)中間點的系數(shù)不能等于0,即因此

3.穩(wěn)定性分析

下面采用von Neumann方法分析本文所推導(dǎo)的差分格式(2.15)的穩(wěn)定性.

對于(2.15)式,舍掉非齊次項f,即假設(shè)f項精確成立,令uni=ηneIσxi,其中,η為振幅,σ為波數(shù),為虛數(shù)單位,有

兩邊同時約掉eIσxi,并整理可得：

利用Euler公式,即eIσh=cosσh+Isinσh,e?Iσh=cosσh ?Isinσh,可得:

對上式進行化簡整理有

從而可得格式(2.15)的誤差放大因子為:

由von Numann穩(wěn)定性定理可知當(dāng)|G| ≤1 時,格式是穩(wěn)定的,由(3.5) 可得|G|=1,因此,格式(2.15) 是無條件穩(wěn)定的.

4.數(shù)值實驗

為了驗證本文格式(2.15)的精確性和穩(wěn)定性,現(xiàn)考慮以下三個具有精確解的初邊值問題.分別采用Crank-Nicolson(C-N)格式,文[7] 中格式和本文格式(2.15) 進行計算; 其中,最大絕對誤差及收斂階的定義為:

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

問題1[7]:

該問題的精確解為:u(x,t)=sin[π(x ?t)].

表1 問題1 當(dāng)λ=τ/h=0.5,t=1 時刻的最大絕對誤差及收斂階

表2 問題1 當(dāng)τ=λh,t=2 時刻的最大絕對誤差

圖1 問題1 當(dāng)N=32,τ=0.03125,t=0.2 時刻的數(shù)值解與精確解

表1給出了針對問題1三種格式在不同空間步長h下,當(dāng)λ=τ/h=0.5,t=1時的最大絕對誤差和收斂階.我們發(fā)現(xiàn)C-N格式在時間和空間上都為二階精度,由于文[7]格式時間具有二階精度,空間具有四階精度,因此當(dāng)取τ=O(h) 時,格式空間僅有二階精度,而本文格式時間和空間均為四階精度.圖1給出N=32,τ=0.03125,t=0.2數(shù)值解與精確解對比圖,可以看出數(shù)值解與精確解吻合的很好.

表2給出了當(dāng)h=1/16和h=1/32時,τ=λh,t=2時刻對問題1采用三種格式計算的最大絕對誤差.可以看出網(wǎng)格比λ最大取到12.8,計算仍然是穩(wěn)定的,因此本文格式是無條件穩(wěn)定的.并且本文格式在所有參數(shù)下,其計算結(jié)果比C-N 格式和文[7]格式計算結(jié)果更加精確.

問題2[7]:

該問題的精確解為:u(x,t)=ecos[π(x?t)].

表3 問題2 當(dāng)λ=τ/h=0.5,t=1 時刻的最大絕對誤差及收斂階

表4 問題2 當(dāng)τ=λh,t=2 時刻的最大絕對誤差

表3和表4給出了針對問題2利用本文格式和C-N格式以及文[7]格式的計算結(jié)果.表3考察了格式的精度,表4驗證了格式的穩(wěn)定性.可以看出本文格式在時間和空間上均可達到四階精度,并且是無條件穩(wěn)定的.

問題3

該問題的精確解為:u(x,t)=cos[π(x ?t)].

表5 問題3 當(dāng)λ=τ/h=0.5,a=0.5,t=1 時刻的最大絕對誤差及收斂階

問題3為非齊次問題,由于文[7]的方程模型為齊次方程,不能計算非齊次問題,因此該問題我們采用本文格式和C-N進行計算和比較,表5給出了兩種格式在不同空間步長h下,當(dāng)t=1時的最大絕對誤差和收斂階.可以看出當(dāng)λ=τ/h=0.5,a=0.5 時,C-N格式在時間和空間上都為二階精度,而本文格式時間和空間均為四階精度.

5.結(jié)論

本文針對一維對流方程提出了一種兩層隱式高精度緊致差分格式,時間和空間均采用泰勒級數(shù)展開法以及截斷誤差余項修正法進行處理,格式截斷誤差為O(τ4+τ2h2+h4),即該格式在時間和空間上均可達到四階精度.并通過von Neumann方法分析得到該格式為無條件穩(wěn)定的.最后通過三個數(shù)值算例驗證了格式的精確性和穩(wěn)定性.通過上述研究,我們可以得出如下結(jié)論:

1.文獻(如[10-11])中的高精度格式往往是時間多層格式,需要另外構(gòu)造啟動步的計算格式,如果采用低精度格式啟動,必然會影響以后時間層的計算精度.而本文格式僅為兩層格式,無須啟動步的計算,時間即可達到四階精度.

2.文獻(如[1,10-11])中的高精度格式空間方向上往往超過三個網(wǎng)格節(jié)點,導(dǎo)致靠近邊界的內(nèi)點計算困難,需要采用特殊處理,而本文格式僅用到三個網(wǎng)格節(jié)點,可以有效避免這一問題.

3.盡管本文格式要求2,這是本文格式的一個缺陷,但是由于本文格式是無條件穩(wěn)定的,從理論上講可以采用任意網(wǎng)格比,因此可以很容易避開2的條件限制,使得這一缺陷并不太影響格式的使用.

5.本文方法可直接推廣到二維和三維問題中去,我們將另文報道.