對流擴散方程的時間間斷時空有限體積元法

肖宇宇 何斯日古楞 楊凱麗

(1.內(nèi)蒙古大學(xué)數(shù)學(xué)科學(xué)學(xué)院,內(nèi)蒙古呼和浩特 010021;2.呼和浩特民族學(xué)院,內(nèi)蒙古呼和浩特 010051)

§1 引言

時間間斷時空有限元方法的研究發(fā)展于20世紀70年代后期.文獻[1]用時間間斷的時空有限元方法求解非線性薛定諤方程時提出了有限元與拉格朗日插值相結(jié)合的技術(shù),即在理論分析時引入由時間剖分單元In上的Radau點確定的拉格朗日插值多項式及相應(yīng)的高斯Radau積分準則,并應(yīng)用有限元技術(shù),插值多項式的基本性質(zhì)及高斯積分準則的高精度性在弱的時空網(wǎng)格限制條件下給出了數(shù)值解的最優(yōu)L∞(L2)模誤差估計.隨后,文獻[2-6]利用這一技術(shù),分別研究了拋物型方程,奇異微分方程,Sobolev方程等發(fā)展型方程的時間間斷時空有限元方法,進行了收斂性分析.文獻[7]針對一維拋物型方程提出了時間間斷的有限體積元格式,并利用有限體積元和有限差分結(jié)合的技巧證明了數(shù)值解的最優(yōu)L∞(L2)模誤差估計.文獻[7]中數(shù)值格式采用分片線性多項式為時空試探函數(shù)空間,只具有二階收斂精度.

本文將時間間斷Galerkin思想與基于等距節(jié)點下三次Lagrange插值的超收斂有限體積元方法[8]相結(jié)合,以三次Lagrange插值導(dǎo)數(shù)超收斂點為對偶剖分節(jié)點,引入插值投影算子,構(gòu)造一種高精度時空有限體積元格式求解對流擴散問題

其中已知函數(shù)p,q,r,f充分光滑,且p(x)≥pmin＞0,r-q′ ≥0.

相比于文獻[1-6]的方法,本文所構(gòu)造的時空有限體積元格式不僅繼承了有限體積元方法[8-9]的高精度,計算簡單等特點,而且能保持物理量的局部守恒性.相對于文獻[7]的方法,這種方法將時間變量和空間變量統(tǒng)一考慮,在時間和空間兩個方向分別發(fā)揮有限元方法和有限體積元方法的優(yōu)勢,從而達到時間和空間變量的高精度,并且在時間剖分節(jié)點處具有超收斂性.

本文結(jié)構(gòu)如下:§1簡單介紹了時空有限元法的發(fā)展歷史及對流擴散問題;§2構(gòu)造了一種高精度時空有限體積元格式;§3分析了所構(gòu)造格式的最優(yōu)收斂性;§4給出數(shù)值算例驗證了所提格式的可行性和有效性.

本文使用標(biāo)準Sobolev空間Hm([a,b]),范數(shù)||v||m,半范數(shù)|v|m以及L2([a,b])空間及其相應(yīng)的內(nèi)積(v,w).本文中所有的c,ci都是與時間和空間步長無關(guān)的正常數(shù),并且不同地方可能取值不相同.

§2 時間間斷時空有限體積元格式

為構(gòu)建高精度時空體積元格式,首先分析三次Lagrange插值多項式的導(dǎo)數(shù)超收斂節(jié)點.在[-h,h]上,用等距節(jié)點(-h,u(-h)),,(h,u(h))對未知函數(shù)u作Lagrange插值,插值函數(shù)記為Πhu,令,則

將u(-h),,u(h)均在x0處作Taylor展開,可得

建立試探函數(shù)空間.在時空片In×Tnh內(nèi),設(shè)Snh ?H10([a,b])是分片三次Lagrange 有限元空間,即

下面給出收斂性分析所需一些定義和引理.引入Radau積分公式[1]

使得對所有次數(shù)不高于2q-2次的多項式精確成立.

使用線性變換t=tn+sΔtn,可將區(qū)間映射到[0,1],從而有如下記號

設(shè)M,N是q×q矩陣,使得

顯然,M和N不依賴于Δtn,且當(dāng)YT=(yn,1,yn,2……·yn,q)∈Rq時,有

引理1[1]設(shè)矩陣,D=diag{s1,s2,……,sq},則有

定義2時空范數(shù)定義為

引理2(參見[8,9]) 對于任意的uh,vh ∈Snh,有

引理3(參見[8]) 對充分小的h,存在正常數(shù)c3,c4和c5,λ ＜∞,使得

其中常數(shù)β滿足rmin ＞β ＞0,,并

證為分析,首先將其分解成

由于v ∈Sh是三次多項式,并用分部求和公式及不等式-ab ≥-(a2+b2)/2,得

定義3在時間單元In=(tn,tn+1]上定義以Radau點tn,j為節(jié)點的Lagrange插值算子:C()→Pq-1(In)滿足

Si基諧振式光學(xué)微腔陀螺的核心技術(shù)指標(biāo)為極限靈敏度,其主要取決于微腔結(jié)構(gòu)的質(zhì)量均勻性、表面粗糙度以及微腔直徑(D)等結(jié)構(gòu)參數(shù)。目前Si基微腔結(jié)構(gòu)在微米級時,表面粗糙度能達到1 nm以下,已接近材料表面粗糙度的極限,Q值到107左右[8-9],此時,微腔陀螺的極限靈敏度就主要取決于微腔直徑D值[10]。

進一步,定義W:[0,T]([a,b])使得

為簡單起見,仍用W表示W(wǎng)|In.

引理4(參見[1]) 由插值逼近性質(zhì),函數(shù)u和函數(shù)W(x,t)的誤差具有如下估計

§3 收斂性分析及誤差估計

定理1設(shè)u和U分別為(1)和(3)的解,則滿足

其中Nc表示(j=1,……,N -1)的總數(shù).

證利用前面定義的W,將誤差項分裂為e=U -u=(U -W)+(W -u)=θ+ρ.對ρ有引理5 的估計結(jié)果,因此下面只需對θ進行估計.首先,由時間間斷時空有限體積元格式(3) 得基本誤差方程

代入(6)式得

其中記號δq,i=ln,i(tn+1)=ln,i(tn,q)=1 (i=q),否則δq,i=0,(i/=q),

在(7)式中取φ=θn,i,然后關(guān)于i從1到q求和,并用引理3可得

對于上式中的Σ1,Σ2,Σ3有如下估計[1]

由于?x ∈[a,b],是關(guān)于時間變量t的2q-2次多項式,因此

于是用Cauchy-Schwarz不等式得

其中,Mn是與n無關(guān)的常數(shù),其具體取值參見后面步驟.

因此利用引理3,有

把上式代入(11)式的右端,再對n進行迭代,經(jīng)整理得

現(xiàn)對于固定的n,取Mm=M=Nc(n)(m=1,……,n).這里Nc(n)表示的總數(shù),并且當(dāng)Nc(n)=0或1時,取M=2.于是,當(dāng).因此,有不等式

其中M=Nc(n).于是利用(13)式和(4),對n=0,……,N -1,有

設(shè)?m(s)表示單元E=(-1,1)上的m次Legendre多項式.進一步,在E上定義多項式

則當(dāng)m ≥0時φm+1(s)在E內(nèi)有m+1個實根,sα,α=0,1,……,m.特別地sm=1是根.

設(shè)是L2([a,b])投影算子,即對任意的u ∈L2([a,b])有

則當(dāng)u ∈H4([a,b])([a,b])時,滿足如下估計式[11-12]:

定理2設(shè)u是問題(1)的解,U是格式(3)的解,并且初始值U0=u0,則存在不依賴Δt和h的正常數(shù)c使得對0≤n ≤N,

由Legendre正交展開知

進一步,構(gòu)造q-1次多項式W(x,s)使得W(x,1)=ω(x,1),且其余項為

其中是待定系數(shù),而bj是(5) 中的已給系數(shù).顯然R(1)=0.令ρ=u-W,用分部積分,余項R和投影的定義,有

現(xiàn)在要求滿足q-1階線性方程組

它的系數(shù)矩陣有如下結(jié)構(gòu)

現(xiàn)令θ=W -U,又由于問題(1) 的解u滿足(14),有

首先,取φ=θ,積分并用Young不等式和引理2后有

假設(shè)h=O(Δt),則對所有時間單元求和并用有限元空間的逆性質(zhì)和Gronwall引理,得

其次,取φ=(t-tn)θt,則φn+=0并有限元空間的逆性質(zhì)和引理2有

利用局部逆性質(zhì),上式變?yōu)?/p>

于是結(jié)合(18)式,(19)式和(20)式,得

上式利用Gronwall 不等式,再求和后代入(18)式有

最后利用L2(Ω)投影的性質(zhì)和余項R的定義,可得對0≤n ≤N,

§4 數(shù)值算例

為驗證本文所提格式的可行性以及理論分析結(jié)果的合理性,考慮一維對流擴散方程[13]

其精確解u(x,t)=exp(5x-(0.25+0.01π2)t)sin(πx).在數(shù)值計算過程中,試探函數(shù)空間取時間線性,空間三次Lagrange插值為基函數(shù),檢驗函數(shù)空間取時間線性,空間分片常數(shù)為基函數(shù).

表1.1分別給出固定時間步長Δt=0.001,而空間步長分別取h=1/2,1/4,1/8,1/16時,時空誤差||u-U||L∞(L2)和時間節(jié)點誤差‖uN -UN‖L2的計算值和收斂階.從表中數(shù)據(jù)可以看出,當(dāng)步長h折半遞減時,兩種誤差數(shù)據(jù)的收斂階接近四階,吻合理論分析結(jié)果.

表1.1 當(dāng)固定時間剖分Δt=0.001時,空間方向的誤差及收斂階

表1.2中分別給出固定空間步長h=0.001時,而時間步長分別取k=1/5,1/10,1/20,1/40時,時空誤差||u-U||L∞(L2)和時間節(jié)點誤差‖uN -UN‖L2的計算值和收斂階.從表中數(shù)據(jù)可以看出,當(dāng)步長Δt折半遞減時,時空誤差數(shù)據(jù)的收斂階接近二階最優(yōu)收斂,而時間節(jié)點處誤差數(shù)據(jù)的收斂階接近三階,呈現(xiàn)超收斂性,符合理論分析結(jié)果.

表1.2 當(dāng)固定空間剖分h=0.001時,時間方向的誤差及收斂階

表1.3分別給出了固定時間步長Δt=0.002,而空間步長分別取h=1/20,1/40,1/80,1/160時,t=tN時刻對偶剖分節(jié)點處空間導(dǎo)數(shù)模誤差和原始剖分節(jié)點處空間導(dǎo)數(shù)Eorigi-nodes:=模誤差.從實驗數(shù)據(jù)發(fā)現(xiàn)原始剖分節(jié)點處導(dǎo)數(shù)項的收斂階為三階,而對偶剖分節(jié)點,即最佳應(yīng)力節(jié)點處收斂階為四階,具有超收斂性.這些實驗數(shù)據(jù)說明時間間斷時空體積元格式能夠有效求解對流擴散問題,且實驗數(shù)據(jù)與理論分析結(jié)果吻合.

表1.3 當(dāng)固定時間剖分Δt=0.002時,導(dǎo)數(shù)在空間方向的誤差及收斂階

§5 結(jié)論

本文針對一維對流擴散問題提出了時間間斷時空有限體積元格式,利用有限差分和有限體積元法相結(jié)合的技巧證明了格式的L∞(L2)-模最優(yōu)收斂誤差估計,用單元正交分解法證明了格式在時間節(jié)點處的超收斂估計.最后為了說明文中所提格式的可行性,給出了數(shù)值算例.從理論分析和數(shù)值實驗數(shù)據(jù)可知,本文所提出的時間間斷時空有限體積元方法在時間剖分節(jié)點處具有超收斂性,高于文獻[8]中有限體積元法的收斂階.本文所提出的時空方法有別于此前的時間間斷時空有限元法(參見文獻[1-7]),文獻[1-6]所涉及方法是時間間斷而空間連續(xù)的時空格式,文獻[7]的方法是時間間斷有限元,空間一次有限體積元的時空格式,而本文所提出的方法時間間斷有限元,空間三次有限體積元的時空格式,具有計算精度高且能保持物理量的局部守恒性等優(yōu)點,并通過數(shù)值實驗發(fā)現(xiàn)在最佳應(yīng)力點處導(dǎo)數(shù)的最大絕對誤差收斂階高于其他節(jié)點處的收斂階.為此,后續(xù)工作中將分析導(dǎo)數(shù)超收斂估計.