一種非等距網格差分格式求解含雙邊界層的反應擴散問題

杜金月，王彩華，張文逸

（天津師范大學 數學科學學院，天津300387）

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

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

1 二階奇異擾動反應擴散問題

考慮反應擴散兩點邊值問題

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

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

考慮問題（1）～（2）的解的漸近展開式

其中

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

邊界條件為

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

其中

將區間[a，b]等分成m份（m為整數且為4的倍數），節點記為 a=x0，x1，x2，…，xm=b，網格步長h=（b-a）/m.記 u（x）i為微分方程（1）～（2）在節點的精確值，Ui為設計的數值方法得到的近似解，c（x）i、（fx）i簡記為 ci、fi.

在節點x=xi處改寫方程（1）

其中g（xi）=（fxi）-c（xi）u（xi）.因為

將式（15）和式（16）代入式（14），并去掉二階截斷誤差，整理得

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

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

令

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

σ2為常數.

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

為簡潔，將上式改寫為

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

對于等距節點的粗剖分，利用差分格式FD2可得到一組數值預測解.然后選擇其中2個內節點a+，將區間剖分為3個區域：左邊界層區域，解平穩的中間區域和右邊界層區域，如圖1所示.

圖1 區間剖分示意圖Fig.1 Diagram of interval partition

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

算法1

步驟1等距差分格式FD2.

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

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

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

步驟2.1如圖1將區間[a，b]分成3個區域，分別 為和.

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

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

步驟2.4將中間區域擴展為[p，q]，取

將區間[p，q]剖分為m/2+2份，步長為

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

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

2 數值實驗

例1考慮一般二階反應擴散方程的邊值問題

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

表1給出了當ε=10-6時例1在不同等距網格剖分下差分格式FD1和FD2的最大誤差（Err）和收斂階（r）的比較.由表1可見含擬合因子的差分格式FD2的計算精度明顯優于FD1，且FD2有比較穩定的二階收斂性.

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

在格式FD2中，m取200，計算不同ε下的最大誤差.執行算法1，在步驟1中m取為80，得到一組預測值，步驟2中N取為200，得到非等距情形下FD3的差分解，計算不同ε下的最大誤差.以上結果見表2.由表2可見，在剖分總節點數相同的情況下，FD3的計算精度有所提高，因為在邊界層處加密了節點，所以與等距情形相比，FD3更有利于觀察解在邊界層的變化.

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

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

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

單獨將差分格式FD2與差分格式FD3的誤差曲線進行比較，見圖3.

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

由圖3可見，差分格式FD3在邊界層的誤差明顯小于FD2的誤差，這說明非等距剖分的計算效果優于等距剖分.

例2考慮一般二階反應擴散方程的邊值問題x∈（0，1），邊界條件為 u（0）=0，u（1）=0，該問題的精確解為

表3給出了當ε=10-3時例2在不同等距網格剖分下差分格式FD1和FD2的最大誤差和收斂階的比較，可見FD2的計算效果明顯優于FD1，為比較穩定的二階收斂.

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

計算不同ε下FD2和FD3的最大誤差，結果見表4，其中剖分方法同表2.圖4給出了當ε=10-3時例2采用3種差分格式的誤差曲線.由表4和圖4可見，非等距網格差分格式FD3的計算效果明顯優于另2種格式.

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

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

3 結語

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

本文是在二階差分格式的基礎上引入了擬合因子，接下來可考慮在更高階的差分格式上引入擬合因子，以進一步提高數值解的精度.