帶間斷系數(shù)的奇異攝動對流擴散方程的自適應(yīng)移動網(wǎng)格方法

朱賜雯,劉利斌

(南寧師范大學(xué)廣西應(yīng)用數(shù)學(xué)中心,廣西南寧 530100)

§1 引言

具有非光滑系數(shù)的奇異攝動微分方程廣泛存在于流體力學(xué),量子力學(xué),彈性力學(xué),化學(xué)反應(yīng)器理論,氣體多孔電極理論,氣象學(xué),海洋學(xué)等領(lǐng)域.這類問題的一個顯著特點是高階導(dǎo)數(shù)項包含一個攝動參數(shù)ε,且對流源項和右邊源項在區(qū)域內(nèi)部存在間斷點.當ε →0時,這類問題的解在區(qū)域的邊界和內(nèi)部某個區(qū)域內(nèi)會產(chǎn)生劇烈的振蕩,即存在所謂的邊界層和內(nèi)層現(xiàn)象,參見[1].正因為邊界層或內(nèi)點層的出現(xiàn),在很多情況下很難給出這類問題的精確解.迄今為止,針對具有連續(xù)系數(shù)的奇異攝動微分方程(即存在邊界層),常用的數(shù)值方法主要有層適應(yīng)網(wǎng)格方法(即Shishkin網(wǎng)格和Bakhvalov網(wǎng)格)和自適應(yīng)移動網(wǎng)格算法[1].

在奇異攝動微分方程邊界層問題的數(shù)值方法受到廣泛關(guān)注的同時,一些具有間斷系數(shù)的奇異攝動問題參數(shù)一致的數(shù)值方法也引起了廣大學(xué)者的興趣,并開展了一系列的研究.例如Farrell[2-3]等人分別在Shishkin網(wǎng)格和一類包含間斷點的分段均勻網(wǎng)格上采用混合差分格式解決帶間斷系數(shù)的奇異攝動對流擴散方程,并證明了數(shù)值解幾乎是一階一致收斂的.針對帶間斷系數(shù)的雙參數(shù)奇異攝動二階常微分方程的邊值問題,Prabha[4]等人在內(nèi)層區(qū)域和其他區(qū)域使用了五點差分格式和標準中心差分格式的組合形式,得到了誤差幾乎是二階收斂.岑仲迪[5]采用混合差分方法在Shishkin網(wǎng)格下研究一類帶間斷系數(shù)的奇異攝動對流擴散方程,并證明了該方法在最大范數(shù)上是二階一致收斂的.此外,研究帶間斷系數(shù)的微分方程的數(shù)值方法還有Chandru[6]提出的包含不連續(xù)點的混合單調(diào)差分方法,和文獻[7]中的迎風有限差分方法.考慮帶間斷系數(shù)的單參數(shù)奇異攝動拋物型問題或雙參數(shù)的奇異攝動拋物型問題的過程中,學(xué)者們在均勻網(wǎng)格和特殊網(wǎng)格下分別構(gòu)造了迎風有限差分格式和混合差分格式,并證明了數(shù)值格式的一致收斂性[8-11].

眾所周知,網(wǎng)格的設(shè)計除了特殊網(wǎng)格和層適應(yīng)網(wǎng)格外,還有自適應(yīng)移動網(wǎng)格.自適應(yīng)移動網(wǎng)格方法是一種基于后驗誤差估計的處理方法,該方法的關(guān)鍵是通過控制函數(shù)在問題的整個區(qū)間上構(gòu)造自適應(yīng)移動網(wǎng)格,使得邊界層和內(nèi)層區(qū)域分布的網(wǎng)格點較多[12-13].對于相同的奇異攝動問題,由算法構(gòu)造的自適應(yīng)移動網(wǎng)格比特殊網(wǎng)格更有效逼近精確解.如文獻[13]中的作者對具有延遲項的奇異攝動微分方程設(shè)計了一個自適應(yīng)移動網(wǎng)格算法,在他們的實驗結(jié)果對比中明顯看到自適應(yīng)移動網(wǎng)格下的數(shù)值解和收斂階比Shishkin網(wǎng)格要更精確.其他的奇異攝動問題如奇異攝動對流擴散方程或方程組問題,帶參數(shù)的非線性奇異攝動問題等,同樣可以在他們的實驗結(jié)果對比中體現(xiàn)出自適應(yīng)移動網(wǎng)格的優(yōu)越性(參考文獻[14-19]).但迄今為止,用自適應(yīng)移動網(wǎng)格算法去處理帶間斷系數(shù)的奇異攝動微分方程的文獻并沒有出現(xiàn).因此,本文僅在一維情況下研究了帶間斷系數(shù)的奇異攝動微分方程的自適應(yīng)移動網(wǎng)格算法,方程為

其中0<ε ?1是攝動參數(shù),α1,α2是正的常數(shù),?d=?-∪?+=(0,d)∪(d,1).對流項系數(shù)a(x)和函數(shù)f(x)在x=d ∈(0,1)處間斷,在其他區(qū)域上是足夠光滑的.在這些假設(shè)條件下,問題(1.1)有唯一解[2],且在x=d這一點處會出現(xiàn)振動現(xiàn)象,即有一個內(nèi)層.

所以針對方程(1.1)的內(nèi)層問題,在間斷點與網(wǎng)格剖分的中點重合的網(wǎng)格下給出混合有限差分格式,基于解及其解的導(dǎo)數(shù)性質(zhì),解的正則和奇異分量的導(dǎo)數(shù)估計以及極大值原理,給出了相應(yīng)的誤差分析.再利用控制函數(shù)和網(wǎng)格等分布原理,對間斷點前后部分網(wǎng)格分別構(gòu)造了自適應(yīng)網(wǎng)格生成算法.最后給出數(shù)值算例驗證該算法的有效性.

注1.2在本文中,C是一個與ε和N無關(guān)的常數(shù),且在不同的位置取不同的值.并且對于任意的函數(shù)u,ui表示u(xi).

§2 連續(xù)解的性質(zhì)

本節(jié)研究問題(1.1)的解及其解的導(dǎo)數(shù)的一些性質(zhì),即極大值原理,連續(xù)解的穩(wěn)定性及導(dǎo)數(shù)估計.

引理2.1(見[2]) 設(shè)u ∈C∩C2(?d)是問題(1.1)的精確解,滿足u(0)≥0,u(1)≥0,且Lu(x)≥0,?x ∈?d,[u](d)=0,[u′](d)≥0.則有u(x)≥0,?(x)∈

引理2.2(見[2]) 設(shè)u(x)是問題(1.1)的解.其解滿足下界

為了后面誤差分析的方便,首先將問題(1.1)的解在區(qū)域[0,d]上分解為u=vL+wL,和在區(qū)域[d,1]上分解為u=vR+wR,且解的光滑層分量vL,vR滿足

邊界層部分wL,wR滿足

引理2.3(見[5,p691]) 對整數(shù)0≤k ≤3,(2.1)-(2.5)的解vL,vR,wL,wR滿足以下估計

§3 離散格式和自適應(yīng)網(wǎng)格

首先將區(qū)間[0,d]和[d,1]分別分成N/2個小區(qū)間,即得如下非均勻網(wǎng)格

引理3.1(見[5,p694]) 假設(shè)uN是問題(3.1)的數(shù)值解,那么其解存在穩(wěn)定性估計

構(gòu)造自適應(yīng)網(wǎng)格算法的關(guān)鍵是對于給定的非負函數(shù)M(x),找到一組非均勻網(wǎng)格N={0=x0

則M(x)稱為網(wǎng)格控制函數(shù),(3.2)稱為相應(yīng)的網(wǎng)格等分布原理.

為了構(gòu)造出問題(3.1)的自適應(yīng)網(wǎng)格移動算法,選擇控制函數(shù)

§4 誤差估計

于是,為了推出離散格式的局部截斷誤差估計,先給出如下若干引理.

綜上,引理4.2證明完成.

令TL,i,TR,i分別為離散格式(3.1)在區(qū)間[0,d)和[d,1] 上的局部截斷誤差,根據(jù)截斷誤差的定義,整合引理4.1和引理4.2得

當xi ∈[xk2,1],i ≥k2時,有hR,i ≥Cε,由引理3.3 以及(3.9)得ε ≤ChR,i ≤CN-1,根據(jù)截斷誤差定義,再利用中值定理可得

注4.1眾所周知,若方程的條件數(shù)趨于無窮時,則導(dǎo)致矩陣奇異或者病態(tài).因此,還要注意論文中方程的條件數(shù)是否趨于無窮.在文獻[20]中作者推導(dǎo)出具有Dirichlet邊界條件的方程為非齊次時,其方程的條件數(shù)為Cond eff=O(1/h),而當方程為齊次時，推導(dǎo)出的條件數(shù)為Cond eff=O(1),所以邊界條件會影響著條件數(shù).下面引用[20]中的引理進行說明.

引理4.6(見[20]中引理4.1) 設(shè)(3.1)的離散格式的系數(shù)矩陣為b,那么矩陣b存在上界為

引理4.7(見[20]中定理4.2) 在引理4.6的條件下,條件數(shù)Cond_eff滿足

由上述的引理可以看出,當ε很小時,條件數(shù)Cond_eff→∞.但如果在實際計算過程中,通過變化把邊界條件轉(zhuǎn)變?yōu)?時,由引理4.7 得到條件數(shù)Cond_eff≤C‖f‖∞.則條件數(shù)與ε無關(guān),從而不會影響解的穩(wěn)定性.

§5 自適應(yīng)網(wǎng)格生成算法

其中MN=MN-1.

§6 數(shù)值實驗與結(jié)果分析

為了驗證本文自適應(yīng)網(wǎng)格算法的有效性,考慮如下奇異攝動對流擴散方程

由于該問題的精確解沒有給出,采用下面公式來計算數(shù)值解的絕對誤差

當ε=2-3k,k=0,1,···,6,N=2j,j=6,7,···,11時,表1 列出了自適應(yīng)網(wǎng)格算法計算得到的最大誤差和迭代次數(shù),其中表1中每一個ε對應(yīng)的第一行表示最大誤差eN,第二行表示相應(yīng)的收斂階rN,第三行為上述網(wǎng)格生成算法的迭代次數(shù).顯然,由表1 的數(shù)據(jù)可以看出隨著N的逐漸增大,本文提出的自適應(yīng)網(wǎng)格算法能到達一階收斂.為了進一步說明本文自適應(yīng)移動網(wǎng)格算法的優(yōu)勢,將混合差分格式(3.1)分別在Shishkin網(wǎng)格和本文設(shè)計的自適應(yīng)移動網(wǎng)格下計算出誤差后進行比較,在表2 中可以看到兩個網(wǎng)格的最大誤差和收斂階,通過對比明顯看出在相同的ε和N的情況下,本文自適應(yīng)移動網(wǎng)格算法的誤差明顯小于Shishkin網(wǎng)格的誤差.

表1 自適應(yīng)網(wǎng)格計算得到的數(shù)值結(jié)果

除此之外,當N=64,ε取不同值時,圖1為上述自適應(yīng)網(wǎng)格算法的數(shù)值解曲線圖.當N=64,ε=2-8時,圖2為上述自適應(yīng)網(wǎng)格算法數(shù)值解的迭代過程.很明顯由兩個圖可以看出,例子的解在x=d處出現(xiàn)了內(nèi)層,自適應(yīng)網(wǎng)格算法的網(wǎng)格點也明顯向內(nèi)層x=d處移動.

圖1 數(shù)值解曲線圖(N=64)

圖2 網(wǎng)格迭代過程(ε=2-8,N=64)

§7 結(jié)論