一種非等距網(wǎng)格差分格式求解含雙邊界層的反應(yīng)擴(kuò)散問(wèn)題

杜金月，王彩華，張文逸

（天津師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，天津300387）

奇異擾動(dòng)問(wèn)題是指微分方程中的高階導(dǎo)數(shù)項(xiàng)含有一個(gè)小參數(shù)，奇異擾動(dòng)問(wèn)題廣泛應(yīng)用于流體力學(xué)、化學(xué)反應(yīng)現(xiàn)象、控制論以及電子場(chǎng)等領(lǐng)域[1-3].傳統(tǒng)數(shù)值方法求解此類問(wèn)題常會(huì)因發(fā)生數(shù)值振蕩而失效，因此尋找適應(yīng)小參數(shù)的數(shù)值方法尤為重要.此類問(wèn)題通常有2種處理方法.一種是從數(shù)值格式構(gòu)造的角度出發(fā)，如：文獻(xiàn)[4-5]針對(duì)一維定常對(duì)流擴(kuò)散反應(yīng)方程，基于截?cái)嗾`差余項(xiàng)修正思想，分別構(gòu)造了一種混合型緊致差分格式和一種高精度緊致差分格式；文獻(xiàn)[6]通過(guò)對(duì)微分方程系數(shù)常數(shù)化處理和基于余項(xiàng)修正思想，給出了對(duì)流擴(kuò)散方程的指數(shù)型高精度緊致差分格式；文獻(xiàn)[7]延用文獻(xiàn)[6]的思想，給出了變系數(shù)對(duì)流擴(kuò)散反應(yīng)方程的指數(shù)型差分格式.另一種思路是從網(wǎng)格適應(yīng)剖分角度出發(fā)，如：文獻(xiàn)[8]基于非均勻網(wǎng)格的泰勒展開(kāi)，通過(guò)一階二階導(dǎo)數(shù)的高階中心差分格式，得到對(duì)流擴(kuò)散反應(yīng)方程的非等距的差分格式；文獻(xiàn)[9]給出對(duì)流擴(kuò)散方程的一種非等距差分格式，并結(jié)合了幾種不同的非等距網(wǎng)格剖分方法，雖然該方法在等距情形下僅有二階精度，但數(shù)值實(shí)驗(yàn)表明其求解小邊界層問(wèn)題的效果較好.

文獻(xiàn)[10]針對(duì)對(duì)流反應(yīng)擴(kuò)散方程，在基本差分格式的系數(shù)上添加擬合因子，數(shù)值實(shí)驗(yàn)表明，相比無(wú)擬合因子（擬合因子為1）的格式，結(jié)果精度得到顯著提高.但該方法是在等距網(wǎng)格下Dahlquist差分格式的基礎(chǔ)上引入擬合因子，因此導(dǎo)致數(shù)值解誤差分布不勻，在邊界層附近的誤差相對(duì)較大.本文對(duì)該方法進(jìn)行改進(jìn).首先在較粗的等距節(jié)點(diǎn)上采用含擬合因子的Dahlquist差分格式，初步得到一組數(shù)值估計(jì)解；然后選擇其中2個(gè)內(nèi)節(jié)點(diǎn)（分別靠近左右邊界層）將區(qū)間剖分為3個(gè)區(qū)域，在邊界層區(qū)域采用含有擬合因子的差分格式再次計(jì)算，而在中間解平滑的區(qū)域使用二階中心差分格式.數(shù)值實(shí)驗(yàn)表明本文方法比文獻(xiàn)[10]的方法計(jì)算效果更好，尤其更有利于觀察雙邊界層附近的數(shù)值解.

1 二階奇異擾動(dòng)反應(yīng)擴(kuò)散問(wèn)題

考慮反應(yīng)擴(kuò)散兩點(diǎn)邊值問(wèn)題

擴(kuò)散量u定義在有限區(qū)間I=（a，b）上，邊界條件為

其中c（x）、f（x）滿足一定的光滑性條件使問(wèn)題的解存在唯一.為方便，設(shè)c（x）≥c＞0，這個(gè)條件限定了問(wèn)題的邊界層出現(xiàn)在端點(diǎn)兩側(cè).

考慮問(wèn)題（1）～（2）的解的漸近展開(kāi)式

其中

式（5）為關(guān)于問(wèn)題（1）的退化解（即小參數(shù)ε為0時(shí)的解）.稱v0（x）為左邊界層校正解，w0（x）為右邊界層校正解，設(shè)其滿足下列常微分方程

邊界條件為

求解方程（6）～（10）， 得v（0τ）和w（0η），代入式（4）得

其中

將區(qū)間[a，b]等分成m份（m為整數(shù)且為4的倍數(shù)），節(jié)點(diǎn)記為 a=x0，x1，x2，…，xm=b，網(wǎng)格步長(zhǎng)h=（b-a）/m.記 u（x）i為微分方程（1）～（2）在節(jié)點(diǎn)的精確值，Ui為設(shè)計(jì)的數(shù)值方法得到的近似解，c（x）i、（fx）i簡(jiǎn)記為 ci、fi.

在節(jié)點(diǎn)x=xi處改寫方程（1）

其中g(shù)（xi）=（fxi）-c（xi）u（xi）.因?yàn)?/p>

將式（15）和式（16）代入式（14），并去掉二階截?cái)嗾`差，整理得

記式（17）為差分格式FD1.

為確定擬合因子σ1，考慮式（1）右端為齊次的情況，即令（fx）≡0，則式（4）中u（0x）=0，且令右邊界層校正解w（0η）=0，則有

令

其中A由式（12）確定.分別將 Ui-1、 Ui、 Ui+1代入式（18），因右端齊次時(shí)有式（18）中fi≡0，且左邊界點(diǎn)處的函數(shù)值已知，考慮邊界層附近h→0，則可得到

σ2為常數(shù).

將擬合因子 σk，k=1、2代入式（18），整理得

為簡(jiǎn)潔，將上式改寫為

式（24）～式（26）為帶擬合因子的差分格式，其系數(shù)矩陣為三對(duì)角矩陣，可采用追趕法求解.以下記該等距差分格式為FD2.

對(duì)于等距節(jié)點(diǎn)的粗剖分，利用差分格式FD2可得到一組數(shù)值預(yù)測(cè)解.然后選擇其中2個(gè)內(nèi)節(jié)點(diǎn)a+，將區(qū)間剖分為3個(gè)區(qū)域：左邊界層區(qū)域，解平穩(wěn)的中間區(qū)域和右邊界層區(qū)域，如圖1所示.

圖1 區(qū)間剖分示意圖Fig.1 Diagram of interval partition

在兩邊界層區(qū)域布置更多的網(wǎng)格節(jié)點(diǎn)，而在中間區(qū)域采用較粗的網(wǎng)格.在邊界層區(qū)域仍采用含有擬合因子的差分格式再次計(jì)算，而在中間解平滑的區(qū)域，利用二階中心差分格式.為提高計(jì)算精度，中間區(qū)域可進(jìn)一步擴(kuò)大到（p，q）， 使得（p，q），邊界p、q處的解值可由前面兩邊界層區(qū)域內(nèi)的計(jì)算值確定.記此時(shí)的非等距差分格式為FD3，具體算法如下：

算法1

步驟1等距差分格式FD2.

先將圖1區(qū)間[a，b]等分成m份（m為整數(shù)，且為4的倍數(shù)），， 在區(qū)間應(yīng)用差分格式（24）～（25）求出.在區(qū)間應(yīng)用差分格式（24）和（26）求出.綜上得到一組.

步驟2非等距差分格式FD3.

先記N為區(qū)間[a，b]的非等距剖分的份數(shù)，N＞m，且為偶數(shù).

步驟2.1如圖1將區(qū)間[a，b]分成3個(gè)區(qū)域，分別 為和.

步驟2.2將左區(qū)域細(xì)化剖分為等份，其中步長(zhǎng)為，應(yīng)用差分格式（24）（k=1），得到 Ui，i=0，1，2，…，（N-m/2）/2，其中該區(qū)域左邊界條件為U0=u（a）=α，右邊界條件為.

步驟2.3將右區(qū)域細(xì)化剖分為等份，其中步長(zhǎng)為，應(yīng)用差分格式（24）（k=2），得到 Ui，i=（N+m/2）/2，（N+m/2）/2+1，…，N，其中該區(qū)域左邊界條件為，右邊界條件為 UN=u（b）= β.

步驟2.4將中間區(qū)域擴(kuò)展為[p，q]，取

將區(qū)間[p，q]剖分為m/2+2份，步長(zhǎng)為

應(yīng)用二階中心差分格式，得到Ui，i=（N-m/2）/2-1，（N-m/2）/2，…，（N+m/2）/2+1，其中該區(qū)間的邊界條件為 u（p） =U(N-m/2)/2-1，u（q） =U(N+m/2)/2+1.

綜上得到 Ui，i=0，1，2，…，N.

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

例1考慮一般二階反應(yīng)擴(kuò)散方程的邊值問(wèn)題

x∈（-1，1），邊界條件為 u（-1）=0，u（1）=0，該問(wèn)題的精確解為

表1給出了當(dāng)ε=10-6時(shí)例1在不同等距網(wǎng)格剖分下差分格式FD1和FD2的最大誤差（Err）和收斂階（r）的比較.由表1可見(jiàn)含擬合因子的差分格式FD2的計(jì)算精度明顯優(yōu)于FD1，且FD2有比較穩(wěn)定的二階收斂性.

表1 當(dāng)ε=10-6時(shí)，例1采用FD1和FD2的最大誤差和收斂階Tab.1 Maximum absolute errors and convergence rates of FD1 and FD2 for example 1 when ε=10-6

在格式FD2中，m取200，計(jì)算不同ε下的最大誤差.執(zhí)行算法1，在步驟1中m取為80，得到一組預(yù)測(cè)值，步驟2中N取為200，得到非等距情形下FD3的差分解，計(jì)算不同ε下的最大誤差.以上結(jié)果見(jiàn)表2.由表2可見(jiàn)，在剖分總節(jié)點(diǎn)數(shù)相同的情況下，F(xiàn)D3的計(jì)算精度有所提高，因?yàn)樵谶吔鐚犹幖用芰斯?jié)點(diǎn)，所以與等距情形相比，F(xiàn)D3更有利于觀察解在邊界層的變化.

表2 ε取不同值時(shí)，例1采用FD2和FD3的最大誤差Tab.2 Maximum absolute errors of FD2 and FD3 for example 1 with different values of ε

圖2給出了當(dāng)ε=10-6時(shí)例1采用3種差分格式的誤差曲線，可以看到中心差分格式在邊界層誤差比較大.

圖2 當(dāng)ε=10-6時(shí)，例1采用3種差分格式的誤差曲線Fig.2 Curves of errors of three difference schemes for example 1 when ε=10-6

單獨(dú)將差分格式FD2與差分格式FD3的誤差曲線進(jìn)行比較，見(jiàn)圖3.

圖3 當(dāng)ε=10-6時(shí)，例1采用FD2和FD3的誤差曲線Fig.3 Curves of errors of FD2 and FD3 for example 1 when ε=10-6

由圖3可見(jiàn)，差分格式FD3在邊界層的誤差明顯小于FD2的誤差，這說(shuō)明非等距剖分的計(jì)算效果優(yōu)于等距剖分.

例2考慮一般二階反應(yīng)擴(kuò)散方程的邊值問(wèn)題x∈（0，1），邊界條件為 u（0）=0，u（1）=0，該問(wèn)題的精確解為

表3給出了當(dāng)ε=10-3時(shí)例2在不同等距網(wǎng)格剖分下差分格式FD1和FD2的最大誤差和收斂階的比較，可見(jiàn)FD2的計(jì)算效果明顯優(yōu)于FD1，為比較穩(wěn)定的二階收斂.

表3 當(dāng)ε=10-3時(shí)，例2采用FD1和FD2的最大誤差和收斂階Tab.3 Maximum absolute errors and convergence rates of FD1 and FD2 for example 2 when ε=10-3

計(jì)算不同ε下FD2和FD3的最大誤差，結(jié)果見(jiàn)表4，其中剖分方法同表2.圖4給出了當(dāng)ε=10-3時(shí)例2采用3種差分格式的誤差曲線.由表4和圖4可見(jiàn)，非等距網(wǎng)格差分格式FD3的計(jì)算效果明顯優(yōu)于另2種格式.

表4 ε取不同值時(shí)，例2采用FD2和FD3的最大誤差Tab.4 Maximum absolute errors of FD2 and FD3 for example 2 with different values of ε

圖4 當(dāng)ε=10-3時(shí)，例2采用3種差分格式的誤差曲線Fig.4 Curves of errors of three difference schemes for example 2 when ε=10-3

3 結(jié)語(yǔ)

本文針對(duì)含雙邊界層的微分方程問(wèn)題，將區(qū)間剖分為3部分，左右邊界層采用含擬合因子的差分格式，給出了一種非等距差分格式.由數(shù)值實(shí)驗(yàn)結(jié)果可以看出，在相同份數(shù)的剖分下，非等距差分方法FD3的數(shù)值結(jié)果優(yōu)于等距差分方法FD2，由誤差曲線也可以明顯看出FD3在邊界層附近的誤差較FD2明顯減小.

本文是在二階差分格式的基礎(chǔ)上引入了擬合因子，接下來(lái)可考慮在更高階的差分格式上引入擬合因子，以進(jìn)一步提高數(shù)值解的精度.