層適應(yīng)網(wǎng)格上求解奇異攝動(dòng)問(wèn)題的粒子群算法

周 琴，程立正

1.湖南涉外經(jīng)濟(jì)學(xué)院 信息與機(jī)電工程學(xué)院，長(zhǎng)沙410205

2.湖南師范大學(xué) 計(jì)算與隨機(jī)數(shù)學(xué)教育部重點(diǎn)實(shí)驗(yàn)室，數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，長(zhǎng)沙410081

1 引言

近年來(lái)，在科學(xué)技術(shù)和工程領(lǐng)域中經(jīng)常面臨求解一些帶小參數(shù)的奇異攝動(dòng)問(wèn)題。如為了開采石油、天然氣等資源需要研究地下多孔介質(zhì)中的滲流問(wèn)題；為預(yù)報(bào)天氣狀況需要求解描述大氣運(yùn)動(dòng)的流體力學(xué)和熱力學(xué)方程組等等。奇異攝動(dòng)問(wèn)題描述的現(xiàn)象往往在局部區(qū)域具有奇異性，它的解含有邊界層或內(nèi)層，解或其導(dǎo)數(shù)在此區(qū)域內(nèi)變化非常劇烈。它的解除與變量有關(guān)外，還與攝動(dòng)小參數(shù)有關(guān)。

若在均勻網(wǎng)格上利用數(shù)值方法求解奇異攝動(dòng)問(wèn)題，為了達(dá)到一定的計(jì)算精度，局部的奇異性會(huì)導(dǎo)致求解區(qū)域上的網(wǎng)格過(guò)細(xì)，造成不必要的計(jì)算時(shí)間和數(shù)據(jù)存儲(chǔ)上的浪費(fèi)。并且均勻網(wǎng)格下求解在解的急劇變化區(qū)域會(huì)產(chǎn)生非物理振蕩，得到不滿意的結(jié)果。層適應(yīng)網(wǎng)格[1-3]是能有效求解奇異攝動(dòng)問(wèn)題的一種非均勻網(wǎng)格，包括Shishkin 網(wǎng)格、Bakhvalov 網(wǎng)格、Bakhvalov-Shishkin 網(wǎng)格等，該類網(wǎng)格具有在邊界層局部加密的特性。學(xué)者們?cè)趯舆m應(yīng)網(wǎng)格上研究各類奇異攝動(dòng)問(wèn)題，已取得一些研究成果，如文獻(xiàn)[4-12]。在層適應(yīng)網(wǎng)格的構(gòu)造中，需要選擇網(wǎng)格過(guò)渡點(diǎn)來(lái)確定網(wǎng)格分布。因此，如何選取網(wǎng)格過(guò)渡點(diǎn)中的參數(shù)，使得數(shù)值解能更好地逼近精確解，是非常有意義的一項(xiàng)工作。文獻(xiàn)[13]對(duì)一類奇異攝動(dòng)問(wèn)題，采用差分進(jìn)化算法優(yōu)化參數(shù)進(jìn)行了求解。文獻(xiàn)[14]對(duì)一類含兩個(gè)參數(shù)的奇異攝動(dòng)問(wèn)題，采用差分進(jìn)化算法在Shishkin網(wǎng)格上進(jìn)行了求解。

粒子群算法（PSO）是由Kennedy J和Eberhart R C[15-16]于1995 年共同提出的一種優(yōu)化算法。它從隨機(jī)解出發(fā)，通過(guò)迭代尋找最優(yōu)解，通過(guò)適應(yīng)度來(lái)評(píng)價(jià)解的品質(zhì)。目前粒子群算法已廣泛應(yīng)用于函數(shù)優(yōu)化、神經(jīng)網(wǎng)絡(luò)訓(xùn)練、模糊系統(tǒng)控制等領(lǐng)域。如文獻(xiàn)[17]提出一種基于陀螺儀傳感器結(jié)合改進(jìn)粒子群算法，解決了傳統(tǒng)3D 建模技術(shù)測(cè)量不精確的缺陷。文獻(xiàn)[18]針對(duì)神經(jīng)網(wǎng)絡(luò)權(quán)值選取不精確的問(wèn)題，提出改進(jìn)的粒子群優(yōu)化算法結(jié)合神經(jīng)網(wǎng)絡(luò)動(dòng)態(tài)選取權(quán)值的方法。文獻(xiàn)[19]研究了改進(jìn)的粒子群優(yōu)化算法在無(wú)約束優(yōu)化問(wèn)題中的應(yīng)用。文獻(xiàn)[20]將擬合方法與粒子群算法相結(jié)合，提出一種新的混合型粒子群優(yōu)化算法，用于方程參數(shù)估計(jì)，得到了較好的結(jié)果。

本文研究一類奇異攝動(dòng)對(duì)流擴(kuò)散問(wèn)題：

其中，ε 是一個(gè)極小的正參數(shù)，滿足0 ＜ε ?1。假設(shè)函數(shù)a(x),b(x),f(x)充分光滑，且存在常數(shù)α 使得a(x)≥2α ＞0,b(x)≥0。這類奇異攝動(dòng)問(wèn)題的解在右邊界層變化劇烈。將對(duì)該類奇異攝動(dòng)問(wèn)題構(gòu)造Bakhvalov-Shishkin網(wǎng)格，采用粒子群算法對(duì)網(wǎng)格參數(shù)α 進(jìn)行優(yōu)化，使得所求數(shù)值解的誤差范數(shù)達(dá)到最優(yōu)。

2 Bakhvalov-Shishkin網(wǎng)格上的差分格式

設(shè)網(wǎng)格剖分?jǐn)?shù)N 為偶數(shù)，網(wǎng)格節(jié)點(diǎn)為{xi}(i=0,1,…,N)，定義網(wǎng)格過(guò)渡點(diǎn), Bakhvalov-Shishkin 網(wǎng)格（簡(jiǎn)記為B-S 網(wǎng)格）的構(gòu)造為將求解區(qū)間[0，1]分解為兩個(gè)子區(qū)間[0,1-τ]和[1-τ,1]，再將子區(qū)間[0,1-τ]作等分，將子區(qū)間[1-τ,1]也分為份，但要使得exp(-α(1-x)/ε)等分布在網(wǎng)格節(jié)點(diǎn)xi(i=N/2,N/2+1,…,N)上，得到如下網(wǎng)格分布[12]：

顯然，當(dāng)小參數(shù)ε 和網(wǎng)格剖分?jǐn)?shù)N 確定時(shí)，網(wǎng)格過(guò)渡點(diǎn)τ 由參數(shù)α 決定。選取不同的參數(shù)α 會(huì)使得網(wǎng)格點(diǎn)發(fā)生變化，從而影響求解的精度。

在式（2）定義的B-S網(wǎng)格{xi}(i=0,1,…,N)上，記uNi為xi處對(duì)應(yīng)的數(shù)值解，可構(gòu)造問(wèn)題（1）的迎風(fēng)差分格式如下：

文獻(xiàn)[12]證明了在B-S網(wǎng)格（2）上利用差分格式（3）求解奇異攝動(dòng)問(wèn)題（1）的數(shù)值解一致一階收斂于真解，即：

其中，u(xi)和分別為網(wǎng)格點(diǎn)xi處對(duì)應(yīng)的精確解與數(shù)值解。

3 優(yōu)化網(wǎng)格參數(shù)的粒子群算法

記uN為當(dāng)網(wǎng)格剖分為N 個(gè)子區(qū)間時(shí)在B-S網(wǎng)格（2）上求解奇異攝動(dòng)問(wèn)題（1）的數(shù)值解，它可視為網(wǎng)格參數(shù)α 的函數(shù)uN(α)。由于很多實(shí)際問(wèn)題中很難得到奇異攝動(dòng)問(wèn)題的精確解，因此使用雙重網(wǎng)格的數(shù)值解來(lái)估計(jì)數(shù)值解誤差，記：

其中，u2N(α)為將網(wǎng)格剖分為2N 個(gè)子區(qū)間時(shí)的數(shù)值解，此時(shí)的網(wǎng)格點(diǎn)由網(wǎng)格剖分為N 個(gè)子區(qū)間時(shí)的B-S網(wǎng)格點(diǎn)及相鄰網(wǎng)格的中點(diǎn)組成。

為了使誤差盡可能得小，構(gòu)造如下的目標(biāo)函數(shù)：

以誤差范數(shù)（4）為目標(biāo)函數(shù)尋找最優(yōu)的網(wǎng)格參數(shù)α。

粒子群算法將參數(shù)α 初始化為一群隨機(jī)粒子，通過(guò)迭代找到最優(yōu)解。在每一次迭代中，粒子通過(guò)跟蹤兩個(gè)最優(yōu)解來(lái)更新自己。第一個(gè)是粒子本身所找到的最優(yōu)解，即個(gè)體最優(yōu)解，記為pbest 。另一個(gè)是整個(gè)種群目前找到的最優(yōu)解，即全局最優(yōu)解，記為gbest 。記Vt和xt分別為第t 次迭代前粒子的速度和位置，粒子根據(jù)如下的公式來(lái)更新速度和位置：

其中，ω 為慣性權(quán)重，r1,r2是區(qū)間（0，1）的隨機(jī)數(shù)，c1,c2是學(xué)習(xí)因子。

下面給出采用粒子群算法結(jié)合迎風(fēng)差分格式在Bakhvalov-Shishkin網(wǎng)格上求解奇異攝動(dòng)對(duì)流擴(kuò)散問(wèn)題（1）的算法步驟：

步驟1 對(duì)參數(shù)α 進(jìn)行種群隨機(jī)初始化，由式（2）得到對(duì)應(yīng)的B-S網(wǎng)格。此時(shí)迭代次數(shù)k=0。

步驟2 利用追趕法等對(duì)B-S網(wǎng)格上的差分格式（3）進(jìn)行求解，得到B-S 網(wǎng)格上N+1 個(gè)網(wǎng)格點(diǎn)上的數(shù)值解uN(α)。在原B-S 網(wǎng)格上添加相鄰網(wǎng)格點(diǎn)的中點(diǎn)，得到2N+1個(gè)網(wǎng)格點(diǎn)，再求解差分格式（3），得到數(shù)值解u2N(α)，進(jìn)而計(jì)算出每個(gè)粒子對(duì)應(yīng)的目標(biāo)函數(shù)（4）的適應(yīng)值。

步驟3 對(duì)每個(gè)粒子α，將其當(dāng)前適應(yīng)值與其個(gè)體歷史最佳位置(pbest)對(duì)應(yīng)的適應(yīng)值做比較，如果當(dāng)前的適應(yīng)值更優(yōu)，則將用當(dāng)前位置更新歷史最佳位置pbest。

步驟4 對(duì)每個(gè)粒子α，將其當(dāng)前適應(yīng)值與全局最佳位置(gbest)對(duì)應(yīng)的適應(yīng)值做比較，如果當(dāng)前的適應(yīng)值更高，則將用當(dāng)前粒子的位置更新全局最佳位置gbest。

步驟5 利用式（5）更新粒子α 的速度和位置。由α的新位置，計(jì)算對(duì)應(yīng)的Bakhvalov-Shishkin網(wǎng)格。

步驟6 若迭代次數(shù)k 小于設(shè)定的最大迭代次數(shù)Gmax，轉(zhuǎn)步驟2，并令k=k+1。否則，停止迭代，得到最優(yōu)網(wǎng)格參數(shù)α 及最優(yōu)誤差范數(shù)。

4 數(shù)值實(shí)驗(yàn)與結(jié)果分析

奇異攝動(dòng)對(duì)流擴(kuò)散問(wèn)題（1）中，ε 是擴(kuò)散項(xiàng)系數(shù)，a(x)是對(duì)流項(xiàng)系數(shù)，當(dāng)ε 極小時(shí)，該問(wèn)題為穩(wěn)態(tài)的對(duì)流占優(yōu)問(wèn)題，可以用來(lái)描述穩(wěn)態(tài)的熱傳導(dǎo)、粒子擴(kuò)散、化學(xué)反應(yīng)等現(xiàn)象。下面選擇一個(gè)變系數(shù)算例和一個(gè)常系數(shù)算例進(jìn)行數(shù)值模擬。

4.1 變系數(shù)奇異攝動(dòng)問(wèn)題

考慮奇異攝動(dòng)問(wèn)題：

根據(jù)上章給出的算法步驟，采用粒子群算法在B-S網(wǎng)格上對(duì)該問(wèn)題進(jìn)行數(shù)值實(shí)驗(yàn)。實(shí)驗(yàn)時(shí)取種群規(guī)模為50，最大迭代次數(shù)Gmax為50 次，慣性權(quán)重ω=0.8,學(xué)習(xí)因子c1=c2=1。參數(shù)α 需要滿足的條件為a(x)≥2α ＞0,其中a(x)是問(wèn)題中u'(x)項(xiàng)的系數(shù)，故此例中參數(shù)α 取值范圍應(yīng)滿足。若α 過(guò)小，網(wǎng)格分布（2）中右邊界層集中的網(wǎng)格點(diǎn)越少，不利于反映解的性質(zhì)，故在實(shí)驗(yàn)時(shí)取粒子α 的搜索范圍為[0.05，0.25]。

4.2 常系數(shù)奇異攝動(dòng)問(wèn)題

考慮奇異攝動(dòng)問(wèn)題：

此例中u'(x) 項(xiàng)的系數(shù)為2，網(wǎng)格參數(shù)α 應(yīng)滿足條件0 ＜2α ≤2,實(shí)驗(yàn)時(shí)取參數(shù)α 的搜索范圍為[0.1，1]，其他實(shí)驗(yàn)參數(shù)的選擇同上例。

4.3 實(shí)驗(yàn)結(jié)果及分析

針對(duì)兩個(gè)算例，取網(wǎng)格剖分?jǐn)?shù)N =32，64，128，256，小參數(shù)ε=10-2,10-3,10-4,10-5,10-6進(jìn)行實(shí)驗(yàn)。表1 和表2分別列出了以式（4）為目標(biāo)函數(shù)時(shí)，算例1和算例2中粒子群算法計(jì)算出的最優(yōu)網(wǎng)格參數(shù)及目標(biāo)函數(shù)值?？梢娫谀繕?biāo)函數(shù)最優(yōu)時(shí)，誤差仍保持收斂性。

圖1作出了兩個(gè)算例中N=32 ε=10-4時(shí)，隨著迭代次數(shù)增加目標(biāo)函數(shù)的變化?？梢娊?jīng)過(guò)10次左右的迭代，目標(biāo)函數(shù)可接近最優(yōu)值。針對(duì)兩個(gè)算例的求解，粒子群算法具有較快的收斂速度。

將參數(shù)α 取兩個(gè)固定值（分別為搜索范圍的上限和下限）與α 取粒子最優(yōu)解時(shí)的誤差范數(shù)目標(biāo)函數(shù)（4）作比較，發(fā)現(xiàn)誤差的區(qū)別主要體現(xiàn)在x=1 附近的邊界層。問(wèn)題的解在x=1 附近的邊界層變化劇烈，為了便于觀察邊界層誤差情況，圖2、圖3分別作出了兩個(gè)算例中ε=10-4時(shí)參數(shù)α 取粒子最優(yōu)解與固定值的邊界層誤差曲線對(duì)比，圖（a）和圖（b）分別為網(wǎng)格剖分?jǐn)?shù)N=32 和N=64 的情況。

表1 算例1中粒子群算法計(jì)算出的最優(yōu)網(wǎng)格參數(shù)及目標(biāo)函數(shù)值

表2 算例2中粒子群算法計(jì)算出的最優(yōu)網(wǎng)格參數(shù)及目標(biāo)函數(shù)值

圖1 N=32,ε=10-4 時(shí)的迭代曲線

圖2 算例1中ε=10-4 時(shí)粒子最優(yōu)解與固定值的邊界層誤差曲線對(duì)比

圖3 算例2中ε=10-4 時(shí)粒子最優(yōu)解與固定值的邊界層誤差曲線對(duì)比

表3、表4 分別列出了算例1 中ε=10-2和ε=10-4時(shí)，采用粒子群算法選取網(wǎng)格參數(shù)和取參數(shù)α 為固定值0.05 和0.25 時(shí)的誤差。表5、表6 分別列出了算例2 中ε=10-2和ε=10-4時(shí)，采用粒子群算法選取網(wǎng)格參數(shù)和取參數(shù)α 為固定值0.1和1時(shí)的誤差。

表3 算例1中ε=10-2 時(shí)的誤差比較

誤差曲線和表中數(shù)據(jù)均反映出采用粒子群算法確定的最優(yōu)參數(shù)α 對(duì)應(yīng)的誤差明顯優(yōu)于α 取固定值時(shí)的誤差。特別是與α 取粒子搜索下限相比，采用粒子群算法得到的數(shù)值解精度要高出許多。因此，與人為選擇B-S網(wǎng)格中的網(wǎng)格參數(shù)相比，利用粒子群算法優(yōu)化選擇網(wǎng)格參數(shù)進(jìn)而求解奇異攝動(dòng)問(wèn)題能得到較好的數(shù)值結(jié)果，特別是邊界層的數(shù)值解精度能得到較大提高。

表4 算例1中ε=10-4 時(shí)的誤差比較

表5 算例2中ε=10-2 時(shí)的誤差比較

表6 算例2中ε=10-4 時(shí)的誤差比較

5 結(jié)束語(yǔ)

本文采用粒子群算法優(yōu)化選擇網(wǎng)格參數(shù)，提出了在Bakhvalov-Shishkin網(wǎng)格上求解一類奇異攝動(dòng)問(wèn)題的粒子群算法。實(shí)驗(yàn)結(jié)果表明，與選擇固定的網(wǎng)格參數(shù)相比，采用粒子群算法計(jì)算能得到更好的計(jì)算精度，并且數(shù)值結(jié)果具有收斂性。本文提到的方法也可以推廣到高維問(wèn)題的計(jì)算中。