求解對流擴(kuò)散方程的4階緊致差分格式

李冉冉,王紅玉,開依沙爾·熱合曼

(新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,新疆 烏魯木齊 830046)

0 引言

對流擴(kuò)散方程是在偏微分方程中一個(gè)重要的分支，被用于描述排水膜中的傳熱、溶解物質(zhì)的分散、液體溶質(zhì)的擴(kuò)散等現(xiàn)象，被廣泛應(yīng)用于環(huán)境科學(xué)，流體力學(xué)等眾多領(lǐng)域中.本文研究1維對流擴(kuò)散方程的初邊值問題:

(1)

其中α為擴(kuò)散系數(shù),β為對流系數(shù),f、g1、g2均為已知充分光滑的函數(shù),u為待求未知量.

近年來，求解該問題的數(shù)值方法有有限差分法[1-10]、有限元法[11]和有限體積法[12]等.由于有限差分法十分簡便且有效，因此在求解該問題時(shí)得到很多研究者的青睞.文獻(xiàn)[3]用梯形方法構(gòu)造該問題的時(shí)間2階、空間4階精度的差分格式.文獻(xiàn)[4]對該問題應(yīng)用時(shí)間3階、空間2階精度的差分格式，該格式在數(shù)值模擬上取得較好的結(jié)果，但在空間方向上的精度有待提高.針對這種情況，學(xué)者們對空間變量提出高階緊致差分格式，文獻(xiàn)[5]給出關(guān)于1階、2階導(dǎo)數(shù)4階緊致離散格式并應(yīng)用于兩點(diǎn)邊值問題的數(shù)值計(jì)算.文獻(xiàn)[6-7]將文獻(xiàn)[5]中4階緊致格式應(yīng)用于Kuramoto-Sivashinsky和耦合Burgers方程組的數(shù)值模擬中.文獻(xiàn)[8]對定常對流擴(kuò)散反應(yīng)方程提出4階精度的差分格式，并通過結(jié)合Richardson技術(shù)達(dá)到提高計(jì)算精度的目的.文獻(xiàn)[9]通過結(jié)合緊致差分格式和4階邊值法構(gòu)造出時(shí)間空間同時(shí)具有4階精度的格式.文獻(xiàn)[10]利用3次Hermite插值公式來求解雙曲型電信方程，在時(shí)間方向上取得4階精度.通過以上分析，本文將結(jié)合文獻(xiàn)[5]和文獻(xiàn)[10]的方法，對帶有周期和Dirichlet邊界條件的對流擴(kuò)散方程構(gòu)造時(shí)間空間同時(shí)具有4階精度的絕對穩(wěn)定的差分格式.

本文給出了空間離散格式以及描述了時(shí)間離散算法，并對格式的穩(wěn)定性給出了證明和滿足2種邊界條件下的算例的數(shù)值驗(yàn)證.

1 空間離散

1.1 在周期邊界條件下對流擴(kuò)散方程的差分格式

考慮具有周期邊界的對流擴(kuò)散方程:

(2)

為了逼近在對流擴(kuò)散方程中的空間導(dǎo)數(shù)，將計(jì)算域Ω×[0,T]={(x,t)|0≤x≤L,0≤t≤T}劃分為由點(diǎn)集{(xj,ti)}描述的均勻網(wǎng)格，記xj=jΔx,j=0,1,…,N,Δx=L/N,ti=iΔt,i=0,1,…,M,Δt=T/M,稱Δx和Δt分別為空間步長和時(shí)間步長.

本文的空間變量的1階導(dǎo)數(shù)項(xiàng)和2階導(dǎo)數(shù)項(xiàng)使用4階緊致差分公式[5]

u′j-1+4u′j+u′j+1=3(uj+1-uj-1)/Δx,j=1,2,…,N

(3)

來近似,其中u′j表示u(x)在j處1階導(dǎo)數(shù)的近似值.

式(3)寫成矩陣形式:

L1U′=M1U,

U=(u1,u2,…,uN)T,U′=(u′1,u′2,…,u′N)T.

則空間變量的1階導(dǎo)數(shù)可以表示為

(4)

若u″j表示u(x)在j處的2階導(dǎo)數(shù)的近似值,則有

u″j-1+10u″j+u″j+1=12(uj-1-2uj+uj+1)/Δx2,j=1,2,…,N.

(5)

式(5)寫成矩陣形式:

L2U″=M2U,

U″=(u″1,u″2,…,u″N)T.

則空間變量的2階導(dǎo)數(shù)可以表示為

(6)

將式(4)和(6)代入方程(2)得到常微分方程

dU/dt=A1U,

(7)

1.2 在Dirichlet邊界條件下對流擴(kuò)散方程的差分格式

為了逼近在對流擴(kuò)散方程中的空間導(dǎo)數(shù)，將計(jì)算域Ω×[0,T]={(x,t)|a≤x≤b,0≤t≤T}劃分為由點(diǎn)集{(xi,tj)}描述的均勻網(wǎng)格，記xj=a+(j-1)Δx,j=1,2,…,N,Δx=(b-a)/(N-1),ti=iΔt,i=0,1,…,M,Δt=T/M.

本文的空間變量的1階和2階導(dǎo)數(shù)同樣用上述4階緊致格式[5]來近似,空間變量的1階導(dǎo)數(shù)項(xiàng)用

4u′2+u′3=(-11u1/12-4u2+6u3-4u4/3+u5/4)/Δx,j=1,

u′j-1+4u′j+u′j+1=3(uj+1-uj-1)/Δx,j=2,3,…,N-1,

u′N-2+4u′N-1=(-uN-4/4+4uN-3/3-6uN-2+4uN-1+11uN/12)/Δx,j=N

(8)

來逼近.式(8)寫成矩陣形式:

LxU′=MxU+b1,

Mx=

b1=11(-u1,0,…,0,uN)T/(12Δx),

U=(u2,u3,…,uN-1)T,

U′=(u′2,u′3,…,u′N-1)T.

則空間變量的1階導(dǎo)數(shù)可以表示為

(9)

空間變量的2階導(dǎo)數(shù)用

14u″2-5u″3+4u″4-u″5=12(u1-2u2+u3)/Δx2,

j=1,

u″j-1+10u″j+u″j+1=12(uj-1-2uj+uj+1)/Δx2,j=2,3,…,N-1,

(10)

來逼近.式(10)寫成矩陣形式:

LxxU″=MxxU+b2,

其中

Lxx=

Mxx=

b2=12(u1,0,…,0,uN)T/Δx2,

U″=(u″2,u″3,…,u″N-1)T.

則空間變量的2階導(dǎo)數(shù)可以表示為

(11)

將式(9)和(11)代入方程(1)得常微分方程

dU/dt=A2U+B1b2-C1b1,

圖1為1/4波長換能器和變幅桿模型，1為前蓋板，2為壓電陶瓷堆，3為后蓋板，4、5為階梯型變幅桿。由于超聲加工屬于輕負(fù)載場合，在設(shè)計(jì)夾心式超聲振子時(shí)，可以忽略負(fù)載對共振頻率的影響，按照空載進(jìn)行計(jì)算。當(dāng)系統(tǒng)共振時(shí)，存在某處振動(dòng)位移為零的節(jié)點(diǎn)。該節(jié)點(diǎn)所在平面稱為波節(jié)面，將波節(jié)面AB設(shè)計(jì)在換能器前蓋板上，截面將超聲換能器分為兩部分，根據(jù)一維傳輸線理論可以分別求得這兩部分的頻率方程[9]：

(12)

2 3次Hermite插值法

3次Hermite插值多項(xiàng)式p(x)[10]可以寫成

p(x)=f(x0)h0(x)+f(x1)h1(x)+f′(x0)·H0(x)+f′(x1)H1(x),

(13)

其中

h0(x)=(1+2(x-x0)/(x1-x0))((x-x1)/

(x0-x1))2,

h1(x)=(1+2(x-x0)/(x0-x1))((x-x0)/

(x1-x0))2,

H0(x)=(x-x0)((x-x1)/(x0-x1))2,

H1(x)=(x-x1)((x-x0)/(x1-x0))2.

在如下定理中，給出了3次Hermite插值多項(xiàng)式(13)的誤差.

定理1[13]設(shè)有2個(gè)不同的點(diǎn)x0、x1,且f∈C4[x0,x1].若p(x)是3次Hermite插值多項(xiàng)式，則在區(qū)間[x0,x1]中的每個(gè)x都對應(yīng)于[x0,x1]中的1個(gè)點(diǎn)ξ，有

f(x)-p(x)=f(4)(ξ)(x-x0)2(x-x1)2/24.

(14)

根據(jù)等式(14)和定理1可以得到積分公式

(15)

在周期邊界下，常微分方程(7)應(yīng)用3次Hermite插值(15),有

因此,式(2)的差分格式為

(16)

同理，常微分方程(12)應(yīng)用3次Hermite插值(15),則在Dirichlet邊界下的對流擴(kuò)散方程(1)的差分格式為

(17)

其中b′1=11(-?u1/?t,0,…,0,?uN-2/?t)T/(12Δx),b′2=12(?u1/?t,0,…,0,?uN-2/?t)T/(Δx)2.

3 穩(wěn)定性分析

引理1[14]設(shè)G=(gij)n×n為嚴(yán)格對角占優(yōu)的，則

1)G是可逆的;

2)若G的所有主對角線項(xiàng)都為正，則G的所有特征值都有正實(shí)部;

3)若G是Hermite矩陣,且G的所有主對角線項(xiàng)都為正,則G的所有特征值都為正實(shí)數(shù).

引理2[15]若z的實(shí)部是非正的，則

|(12+6z+z2)/(12-6z+z2)|≤1.

定理2本文差分格式(16)和(17)的截?cái)嗾`差為o(Δt4+Δx4)且無條件穩(wěn)定.

λj=αλll/λmm-βλl/λm.

|(λE)j|≤1(j=1,2,…,n+1).

在本文格式(17)的齊次情況下，結(jié)構(gòu)與格式(16)是一致的，因此，格式(16)、(17)都是無條件穩(wěn)定的.因?yàn)榫o致離散格式(4)、(6)、(9)、(11)在空間方向上為4階，且3次Hermite插值公式在時(shí)間方向上為4階，因此差分格式(16)、(17)截?cái)嗾`差為o((Δt)4+(Δx)4).

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

為了研究本文格式的可靠性，對滿足2種邊界條件下的算例進(jìn)行數(shù)值計(jì)算，并比較了本文格式的L2范數(shù)誤差及空間時(shí)間的收斂階.表1～表4給出了：本文格式(16)、(17)的時(shí)間變量均為Crank-Nicolson，空間變量分別為中心差分格式[1]o((Δx)2+(Δt)2)、4階差分格式[2](HOC)

o((Δx)4+(Δt)2).其中L2范數(shù)誤差和收斂階的定義分別為

r=lg(‖u-UΔx(Δt)‖/‖u-UΔx(Δt)/2‖/lg 2,

其中Uk、uk分別表示點(diǎn)xk處的數(shù)值解和準(zhǔn)確解.

例1考慮具有周期邊界的對流擴(kuò)散方程[16]:

其準(zhǔn)確解為u(x,t)=e-αtsin(x-βt),其中α=β=1,L=2π.

表1和表2針對周期邊界的對流擴(kuò)散方程,當(dāng)T=1、Δt=0.002和T=1、Δx=1/60時(shí),給出了在不同Δx和Δt的L2范數(shù)誤差和收斂階的比較.由表1和表2可以看出:當(dāng)空間變量離散分別為2階中心差分格式、HOC,時(shí)間變量離散為Crank-Nicolson時(shí),本文格式(16)具有更小的誤差,且在時(shí)間和空間方向上都達(dá)到了4階精度,這與理論結(jié)果相符合.

例2考慮具有Dirichlet邊界的對流擴(kuò)散方程[17]:

表1 例1中當(dāng)T=1、Δt=0.002、不同Δx時(shí)的L2范數(shù)誤差和收斂階比較

表2 例1中當(dāng)T=1、Δx=1/60、不同Δt時(shí)的L2誤差和收斂階比較

表3和表4針對Dirichlet邊界的對流擴(kuò)散方程,當(dāng)T=1、Δt=(Δx)2和T=1、Δx=0.012 5時(shí),給出了在不同Δx和Δt時(shí)的L2范數(shù)誤差和收斂階的比較.由表3和表4可以看出:當(dāng)空間變量離散分別為2階中心差分格式、HOC,時(shí)間變量離散為Crank-Nicolson時(shí),本文格式(17)具有更小的誤差,且在時(shí)間和空間方向上都達(dá)到了4階精度,這與理論結(jié)果相符合.

表3 例2中當(dāng)T=1、Δt=(Δx)2、不同Δx時(shí)的L2范數(shù)誤差和收斂階比較

表4 例2中當(dāng)T=1、Δx=0.012 5、不同Δt時(shí)的L2誤差和收斂階比較

5 結(jié)論

本文在周期和Dirichlet邊界條件下對空間變量的離散達(dá)到了4階精度，并結(jié)合3次Hermite插值法提出了數(shù)值求解1維對流擴(kuò)散方程在空間和時(shí)間上都具有4階精度且絕對穩(wěn)定的差分格式.采用本文的方法計(jì)算了2個(gè)數(shù)值算例，從表1和表3可知:本文的格式在空間上達(dá)到了4階的精度，與傳統(tǒng)的中心差分格式和HOC格式相比，體現(xiàn)了本文格式的有效性；從表2和表4可知:本文的格式在時(shí)間上也達(dá)到了4階的精度，與Crank-Nicolson相比，本文的格式更精確.另外，該格式也可以通過利用局部1維化方法推廣到2維和3維對流擴(kuò)散方程問題的數(shù)值計(jì)算中.