奇異攝動內層問題的寬度估計及其數值求解

邵文婷，鄭爍宇

1.上海第二工業大學 文理學部，上海201209

2.同濟大學 數學科學學院，上海200092

1 引言

在物理和工程技術等應用領域，許多從實際問題抽象得到的數理方程中，最高階導數項前通常含有反映介質性質的攝動參數ε。例如，對流擴散[1-2]和反應擴散方程[3-4]、半導體設備建模的漂移擴散方程、具有大Peclet數的熱傳導方程等。問題的解作為ε的函數具有奇性層現象：當ε趨于0時，在很小的區域內（區域寬度與ε有關），函數值變化非常劇烈，然而在遠離奇性的區域，變化相對比較平緩。這類問題稱為奇異攝動問題。

對于奇異攝動問題一般難以得到解析解，而這類跨尺度攝動現象在應用上普遍存在，又使得其迫切需要高效快速的近似解法。目前，隨著計算機技術高度發展，數值解法引起了廣大科學工作者的重視[5-8]，逐漸活躍為主流方法之一。然而，解具有的奇性層效應給有限差分、有限元等經典的數值方法帶來了挑戰，往往得不到理想的計算精度。因此，眾多學者致力于尋找一種有效的奇性層校正方法，它在實現高精度數值計算上發揮著重大作用。在校正問題的研究中大多涉及了奇性層寬度的探究。文獻[9]提出了一種sinh變換，使得離散節點在奇性區域內部進行加密。sinh變換表達式中含有兩個分別用來刻畫奇性位置和寬度的參數。位置參數可以由漸近理論對解進行先驗分析來預判。對于寬度參數，文獻[9]和文獻[10]分別提出了基于Chebyshev-Páde逼近和Fourier-Páde逼近的估計方法。然而，這些方法在各種實際應用問題中的計算效果一般，可操作性不強。文獻[5-6，11]考慮了一類奇異攝動兩點邊值問題，構造了寬度公式d=σε。在此基礎上，文獻[12]考慮了一類二維規則區域上指數型邊界層問題，提出一種具有“邊界消減技術”的區域分裂方法。思路是基于寬度公式對求解區域進行合理的劃分，在每個子區域上采用Chebyshev-Tau高精度方法分別進行離散求解。然而，對于公式中系數σ的選取，這些工作均依賴于人為經驗，并沒有給出一個具體的計算策略。近年來，劉利斌等[4，13-15]在寬度的估計上進行了深入研究。文獻[4]針對Shinkin網格在求解反應擴散方程的應用上，構造了以誤差范數最小為目標函數的無約束優化問題，采用粒子群優化算法給出了網格過渡點參數的估計。文獻[13]將這類寬度估計方法推廣至含有兩個攝動參數的問題，借助差分進化算法使得邊界層區域的計算精度得到大幅度提高。文獻[14-15]將利用非線性優化問題來估計寬度參數的想法應用于帶有sinh變換的有理譜配點法，加強了sinh變換在求解邊界層問題[14]和內層問題[15]中的可靠性。從這些工作中，不難發現尋求一種切實合理的寬度估計方法可以使得奇異攝動問題的數值求解更加精確、高效。

受上述研究啟發，本文針對奇異攝動內層問題提出了一種內層寬度的估計方法，并且采用有理譜配點法（Rational Spectral Collocation Method，RSCM）進行高精度數值求解。首先，將計算區域分解為內層區域和區域外部，原問題相應地分解為內層問題和退化問題兩部分，其中退化問題不再具有奇性。對于內層問題，引入一個適當的拉伸變換來減弱奇性。為了確定內層寬度，基于區域連接端點處的一階導數連續性條件設計了一個非線性優化問題。在數值計算上，采用重心形式有理譜配點法進行離散。

重心形式有理譜配點法最早是由Berrut[16-17]等人提出的，其構造特點是將函數在各個節點處的導數值由計算區域上所有離散節點處的函數值加權之和來近似，將微分方程轉化為代數方程組進行求解。它的主要優點是如果依據問題的需要引入坐標變換，原方程并不需要隨之改變，易于數值離散的操作實現[17]。在計算精度方面，問題的解是充分光滑的，可以達到指數階收斂速度，因而備受計算數學研究者的關注[5-6，14-15，18-19]。

本文的主要貢獻是基于RSCM構造了一種內層寬度估計的數值方法。該方法結合了區域分解和拉伸變換，設計了關于一階導數連續性條件的優化問題來確定內層寬度，克服了依賴人為經驗設定參數的不合理性，或是采用試錯法和枚舉法造成的計算量和計算時間的大量浪費。此外，本文還考慮了一類非線性奇異攝動兩點邊值問題的數值計算，采用Adomian分解法[20-21]對原問題進行線性化處理，從而可以將已有的研究工作推廣至非線性問題的求解上。數值實驗表明，對線性和非線性問題的內層寬度，本文提出的數值方法均給出了一個合理有效的估計，同時在數值解的計算上也得到了令人滿意的精度。

2 有理譜配點法

其中，lj(x)為插值基函數：

相比經典的拉格朗日型插值函數，重心形式有理插值函數的顯著優點是：當離散節點個數較多或需要求解的線性方程組規模較大時，計算引入的舍入誤差要小得多，數值穩定性更強。

另一方面，通過簡單推導可以直接得到未知函數導數的逼近。設Iu(x)在節點xj處的r階導數為：

其中一階和二階微分矩陣的元素表達式為[9]：

本文取Chebyshev-Gauss-Lobatto點為離散節點，即xj=cos(jπ/N)，對應的權系數為ωj=(-1)jδj,δ0=δN=0.5,δj=1,j=2,…,N-1。

3 區域分解及內層寬度的數值估計

3.1 區域分解

本文考慮如下一類奇異攝動兩點邊值問題：

其中，0＜ε?1為攝動參數，ua和ub為已知邊界條件。f(x,u)和g(x,u)關于x、u一階可微。依據f(x,u)和g(x,u)的取值，問題（6）的解在區間端點出現邊界層或者區間內部出現內層[14-15]。本文僅考慮內層問題的情況。

假設問題的解在點x?∈(a,b)的一個小鄰域出現內層，寬度d與參數ε有關。不妨假設d=2σεm，其中m＞0,σ＞0為待定參數，這里稱σ為寬度系數。

首先，借鑒量綱分析的方法確定m。對問題（6）引入拉伸變換x=εmξ，控制方程轉化為：

為了使得變換后的方程對于新尺度ξ的量綱是一致的，選取m=m?使得下式成立：

隨后，將求解區域分解為兩部分：內層區域[x?-σεm?,x?+σεm?]和外部區域[a,x?-σεm?]，[x?+σεm?,b]。對應原問題（6）的解作如下分解。

問題1在外部區域，滿足退化問題：

以及

退化問題不含有攝動參數ε，是無奇性問題。一旦確定寬度系數σ，便可以采用第2章介紹的RSCM進行離散求解。

問題2在內層區域，滿足內層問題：

其中，uL(x)和uR(x)分別記為退化問題（9）和（10）的解。

對內層問題（11）引入拉伸變換x=x?+σεm?η，η∈[-1,1]，并且定義下列函數：

結合式（8），內層問題（11）轉化為如下邊值問題：

令ε→0，得到內層問題的近似：

隨后，將通過求解該近似問題來確定寬度系數σ。

3.2 寬度估計及離散格式

本節采用RSCM對問題（13）進行離散，同時給出寬度系數σ的一個數值估計。

首先，建立一階和二階微分矩陣的分塊形式：

由矩陣和向量的運算性質，得到問題（13）的離散形式：

其中算子“°”表示矩陣的哈達瑪（Hadamard）積。

需要指出的是，H?是一個不僅與未知量u?2有關，還與寬度系數σ有關的系數矩陣。因此，方程組（15）所含未知量的個數大于方程的個數，是一個不定問題。

為了確定寬度系數σ，同時求得內層問題的數值解，本文基于內層區域和外部區域連接端點處的一階導數連續性條件，構造如下優化問題：

優化問題（16）的離散形式與代數方程組（15）構成一個以u?2與σ為未知量的非線性方程組，這里通過Matlab優化工具箱中非線性方程求解器fsolve來進行計算，其本質原理是最小二乘法。在實際操作中，迭代初始值可以取為邊界條件的線性插值。

將求得的最優寬度系數σ?與式（8）確定的m?代入寬度公式d=，至此便得到了內層寬度的估計值。為了進一步提高計算精度，可以將m?及σ?代入退化問題（9）、（10）和內層問題（12），采用RSCM進行離散求解，得到更高精度的數值解。

4 一類非線性奇異攝動問題的線性化求解

對于非線性問題，通常需要選取合適的迭代法進行求解。當方程的非線性項比較復雜時，牛頓迭代等一些經典格式的求解效果往往不理想，更不用提問題本身還具有奇性。本文引入文獻[20-21]提出的Adomian分解法將非線性方程進行線性化。由數值實驗可以得到驗證，該線性化方法對奇異攝動問題的求解是行之有效的。

以如下一類非線性兩點邊值問題為模型介紹Adomian分解法：

記un-1(n≥2)為第n-1層迭代得到的數值解。在第n層上，對未知函數un進行分解：

同時，導函數有如下分解：

將式（18）～（20）代入模型問題（17），整理得到v0,1,2,3分別滿足下列問題：

注意到un-1是已知的，因此方程（21）～（24）是線性問題，可以直接采用RSCM離散求得v0,1,2,3，代入式（18）得到un。當相鄰兩層迭代得到的數值解滿足給定的容忍誤差＜Tol時，計算終止。

現將Adomian分解線性化迭代算法的求解過程做如下整理：

步驟2當n≥1時，采用RSCM依次離散求解線性方程（21）～（24）得到v0,1,2,3，其中各方程所涉及的導數項可以由式（5）計算得到。將v0,1,2,3代入式（18）得到un。

5 數值實驗

本文給出兩個算例，來驗證數值估計方法的可行性和有效性。算法通過Matlab軟件實現，計算機的運行環境為Intel?CoreTMi5-6500 CPU 3.20 GHz，4.0 GB，Windows7（x64）。

例1考慮線性問題[22]：

具有精確解：

其中，erf(x)為剩余誤差函數，δ=1/2ε。

根據奇異攝動漸近理論，該問題的解在x=0的小鄰域內出現奇性層，不妨假設內層區域為[-σεm,σεm]。

首先，借鑒量綱分析方法確定m?。對控制方程引入拉伸變換x=εmξ，轉化為：

取1-2m=0，得到m?=1/2，即內層寬度為d=2σε。

隨后，將原問題分解為[-1,-σε]和[σε,1]上的退化問題，以及[-σε,σε]上的內層問題。采用第3.2節提出的寬度估計方法與RSCM結合進行離散求解。

取ε=10-3,10-4,…,10-9進行計算，每個區間上離散節點數N=30。由優化問題（16），求得最優寬度系數σ?=4.5。表1通過枚舉法列出了對應不同的ε值，當1≤σ≤7時，RSCM求得的離散節點上的最大誤差。對比這些數值結果，可以發現當4≤σ≤6時，RSCM對所有的ε值都得到了比較理想的計算精度，并且精度值穩定。而遠離該區間的σ值，計算誤差逐步增大，本文方法求得的σ?=4.5落在該區間內，驗證了其合理性和有效性，避免采用枚舉法或試錯法產生的大量計算時間的浪費。圖1刻畫了ε=10-3,10-4,10-5時RSCM求得的數值解，可以看出解的變化趨勢與估計得到的內層寬度是相互吻合的。

表1 例1中ε、σ取不同值RSCM計算得到的最大誤差

圖1 例1中ε=10-3,10-4,10-5時的數值解及內層寬度

例2考慮如下非線性奇異攝動問題：

具有精確解：

由漸近理論分析得到，該問題的解在x=0的小鄰域內出現奇性層，不妨假設內層區域為[-σεm,σεm]。

首先，借鑒量綱分析方法確定m?。對控制方程引入拉伸變換x=εmξ，轉化為：

取m-1=0，得到m?=1，即內層寬度為d=2σε。原問題相應地分解為[-1,-σε]和[σε,1]上的退化問題，以及[-σε,σε]上的內層問題。

本算例中考查ε=10-3,10-4,10-5時的數值求解，每個區間上取離散節點數N=50。由優化問題（16）求得最優寬度系數σ?=8.2。進一步，結合第4章介紹的Adomian分解法迭代計算求得數值解，其中容忍誤差Tol=10-8。表2通過枚舉法列出了5≤σ≤13時，RSCM計算得到的最大誤差。觀察發現當7≤σ≤10時，RSCM得到的計算精度較優，并且穩定在10-4～10-5，本文方法得到的σ?=8.2落在該區間內。當σ落在區間外部時，誤差逐漸增大。當σ≥13時，甚至出現計算不收斂現象。由此可見，內層寬度的估計是否合理對非線性奇異攝動問題的求解起著至關重要的作用。此外，圖2刻畫了ε=10-3,10-4,10-5時，RSCM求得的數值解及估計出的內層寬度。這些數值結果說明了本文提出的數值方法對非線性問題也是有效的，同時得到了不錯的計算精度。

表2 例2中ε、σ取不同值RSCM計算得到的最大誤差

圖2 例2中ε=10-3,10-4,10-5時的數值解及內層寬度

6 結束語

本文考慮了奇異攝動內層問題的數值計算，提出了一種內層寬度的估計方法。首先，將原問題依據內層區域和區域外部相應地分解為內層問題和退化問題兩部分。隨后，構造了一個非線性優化問題來確定寬度系數。采用重心形式的有理譜配點法進行離散求解。對一類非線性奇異攝動問題，引入Adomian分解法進行線性化迭代計算。數值實驗表明，本文提出的方法對線性和非線性問題都是可行的，切實有效地給出了內層寬度的一個合理估計，數值解也達到了不錯的計算精度。相比依賴人為經驗的試錯法或是枚舉法，避免了大量計算時間的耗費，更是一種經濟的數值方法。值得一提的是，本文的數值方法可以直接應用于邊界層問題的計算。

目前，奇異攝動問題的數值解法還有許多問題值得深入研究。例如，求解區域上存在多個奇性層以及高維問題，這些都需要對本文方法進行不斷改進，以適應求解各種不同類型的奇性問題。本文提出的內層寬度估計方法及實驗結果，能為奇異攝動問題的理論分析及其數值計算的發展提供一臂之力。