奇異攝動內(nèi)層問題的寬度估計及其數(shù)值求解

邵文婷，鄭爍宇

1.上海第二工業(yè)大學(xué) 文理學(xué)部，上海201209

2.同濟大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，上海200092

1 引言

在物理和工程技術(shù)等應(yīng)用領(lǐng)域，許多從實際問題抽象得到的數(shù)理方程中，最高階導(dǎo)數(shù)項前通常含有反映介質(zhì)性質(zhì)的攝動參數(shù)ε。例如，對流擴散[1-2]和反應(yīng)擴散方程[3-4]、半導(dǎo)體設(shè)備建模的漂移擴散方程、具有大Peclet數(shù)的熱傳導(dǎo)方程等。問題的解作為ε的函數(shù)具有奇性層現(xiàn)象：當(dāng)ε趨于0時，在很小的區(qū)域內(nèi)（區(qū)域?qū)挾扰cε有關(guān)），函數(shù)值變化非常劇烈，然而在遠離奇性的區(qū)域，變化相對比較平緩。這類問題稱為奇異攝動問題。

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

受上述研究啟發(fā)，本文針對奇異攝動內(nèi)層問題提出了一種內(nèi)層寬度的估計方法，并且采用有理譜配點法（Rational Spectral Collocation Method，RSCM）進行高精度數(shù)值求解。首先，將計算區(qū)域分解為內(nèi)層區(qū)域和區(qū)域外部，原問題相應(yīng)地分解為內(nèi)層問題和退化問題兩部分，其中退化問題不再具有奇性。對于內(nèi)層問題，引入一個適當(dāng)?shù)睦熳儞Q來減弱奇性。為了確定內(nèi)層寬度，基于區(qū)域連接端點處的一階導(dǎo)數(shù)連續(xù)性條件設(shè)計了一個非線性優(yōu)化問題。在數(shù)值計算上，采用重心形式有理譜配點法進行離散。

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

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

2 有理譜配點法

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

相比經(jīng)典的拉格朗日型插值函數(shù)，重心形式有理插值函數(shù)的顯著優(yōu)點是：當(dāng)離散節(jié)點個數(shù)較多或需要求解的線性方程組規(guī)模較大時，計算引入的舍入誤差要小得多，數(shù)值穩(wěn)定性更強。

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

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

本文取Chebyshev-Gauss-Lobatto點為離散節(jié)點，即xj=cos(jπ/N)，對應(yīng)的權(quán)系數(shù)為ωj=(-1)jδj,δ0=δN=0.5,δj=1,j=2,…,N-1。

3 區(qū)域分解及內(nèi)層寬度的數(shù)值估計

3.1 區(qū)域分解

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

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

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

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

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

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

問題1在外部區(qū)域，滿足退化問題：

以及

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

問題2在內(nèi)層區(qū)域，滿足內(nèi)層問題：

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

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

結(jié)合式（8），內(nèi)層問題（11）轉(zhuǎn)化為如下邊值問題：

令ε→0，得到內(nèi)層問題的近似：

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

3.2 寬度估計及離散格式

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

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

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

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

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

為了確定寬度系數(shù)σ，同時求得內(nèi)層問題的數(shù)值解，本文基于內(nèi)層區(qū)域和外部區(qū)域連接端點處的一階導(dǎo)數(shù)連續(xù)性條件，構(gòu)造如下優(yōu)化問題：

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

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

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

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

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

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

同時，導(dǎo)函數(shù)有如下分解：

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

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

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

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

5 數(shù)值實驗

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

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

具有精確解：

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

根據(jù)奇異攝動漸近理論，該問題的解在x=0的小鄰域內(nèi)出現(xiàn)奇性層，不妨假設(shè)內(nèi)層區(qū)域為[-σεm,σεm]。

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

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

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

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

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

圖1 例1中ε=10-3,10-4,10-5時的數(shù)值解及內(nèi)層寬度

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

具有精確解：

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

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

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

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

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

圖2 例2中ε=10-3,10-4,10-5時的數(shù)值解及內(nèi)層寬度

6 結(jié)束語

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

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