邵文婷，鄭爍宇

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

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

1 引言

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

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

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

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

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

2 有理譜配點(diǎn)法

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

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

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

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

本文取Chebyshev-Gauss-Lobatto點(diǎn)為離散節(jié)點(diǎn)，即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ù)值估計(jì)

3.1 區(qū)域分解

本文考慮如下一類奇異攝動(dòng)兩點(diǎn)邊值問題：

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

假設(shè)問題的解在點(diǎn)x?∈(a,b)的一個(gè)小鄰域出現(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ū)域，滿足退化問題：

以及

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

問題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 寬度估計(jì)及離散格式

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

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

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

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

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

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

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

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

4 一類非線性奇異攝動(dòng)問題的線性化求解

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

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

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

同時(shí)，導(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時(shí)，計(jì)算終止。

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

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

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

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

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

具有精確解：

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

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

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

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

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

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

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

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

例2考慮如下非線性奇異攝動(dòng)問題：

具有精確解：

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

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

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

6 結(jié)束語

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

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