帶間斷系數的奇異攝動對流擴散方程的自適應移動網格方法

朱賜雯,劉利斌

(南寧師范大學廣西應用數學中心,廣西南寧 530100)

§1 引言

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

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

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

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

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

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

§2 連續解的性質

本節研究問題(1.1)的解及其解的導數的一些性質,即極大值原理,連續解的穩定性及導數估計.

引理2.1(見[2]) 設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]) 設u(x)是問題(1.1)的解.其解滿足下界

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

邊界層部分wL,wR滿足

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

§3 離散格式和自適應網格

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

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

構造自適應網格算法的關鍵是對于給定的非負函數M(x),找到一組非均勻網格N={0=x0

則M(x)稱為網格控制函數,(3.2)稱為相應的網格等分布原理.

為了構造出問題(3.1)的自適應網格移動算法,選擇控制函數

§4 誤差估計

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

綜上,引理4.2證明完成.

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

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

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

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

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

由上述的引理可以看出,當ε很小時,條件數Cond_eff→∞.但如果在實際計算過程中,通過變化把邊界條件轉變為0時,由引理4.7 得到條件數Cond_eff≤C‖f‖∞.則條件數與ε無關,從而不會影響解的穩定性.

§5 自適應網格生成算法

其中MN=MN-1.

§6 數值實驗與結果分析

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

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

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

表1 自適應網格計算得到的數值結果

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

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

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

§7 結論