層適應網格上求解奇異攝動問題的粒子群算法

周 琴，程立正

1.湖南涉外經濟學院 信息與機電工程學院，長沙410205

2.湖南師范大學 計算與隨機數學教育部重點實驗室，數學與統計學院，長沙410081

1 引言

近年來，在科學技術和工程領域中經常面臨求解一些帶小參數的奇異攝動問題。如為了開采石油、天然氣等資源需要研究地下多孔介質中的滲流問題；為預報天氣狀況需要求解描述大氣運動的流體力學和熱力學方程組等等。奇異攝動問題描述的現象往往在局部區域具有奇異性，它的解含有邊界層或內層，解或其導數在此區域內變化非常劇烈。它的解除與變量有關外，還與攝動小參數有關。

若在均勻網格上利用數值方法求解奇異攝動問題，為了達到一定的計算精度，局部的奇異性會導致求解區域上的網格過細，造成不必要的計算時間和數據存儲上的浪費。并且均勻網格下求解在解的急劇變化區域會產生非物理振蕩，得到不滿意的結果。層適應網格[1-3]是能有效求解奇異攝動問題的一種非均勻網格，包括Shishkin 網格、Bakhvalov 網格、Bakhvalov-Shishkin 網格等，該類網格具有在邊界層局部加密的特性。學者們在層適應網格上研究各類奇異攝動問題，已取得一些研究成果，如文獻[4-12]。在層適應網格的構造中，需要選擇網格過渡點來確定網格分布。因此，如何選取網格過渡點中的參數，使得數值解能更好地逼近精確解，是非常有意義的一項工作。文獻[13]對一類奇異攝動問題，采用差分進化算法優化參數進行了求解。文獻[14]對一類含兩個參數的奇異攝動問題，采用差分進化算法在Shishkin網格上進行了求解。

粒子群算法（PSO）是由Kennedy J和Eberhart R C[15-16]于1995 年共同提出的一種優化算法。它從隨機解出發，通過迭代尋找最優解，通過適應度來評價解的品質。目前粒子群算法已廣泛應用于函數優化、神經網絡訓練、模糊系統控制等領域。如文獻[17]提出一種基于陀螺儀傳感器結合改進粒子群算法，解決了傳統3D 建模技術測量不精確的缺陷。文獻[18]針對神經網絡權值選取不精確的問題，提出改進的粒子群優化算法結合神經網絡動態選取權值的方法。文獻[19]研究了改進的粒子群優化算法在無約束優化問題中的應用。文獻[20]將擬合方法與粒子群算法相結合，提出一種新的混合型粒子群優化算法，用于方程參數估計，得到了較好的結果。

本文研究一類奇異攝動對流擴散問題：

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

2 Bakhvalov-Shishkin網格上的差分格式

設網格剖分數N 為偶數，網格節點為{xi}(i=0,1,…,N)，定義網格過渡點, Bakhvalov-Shishkin 網格（簡記為B-S 網格）的構造為將求解區間[0，1]分解為兩個子區間[0,1-τ]和[1-τ,1]，再將子區間[0,1-τ]作等分，將子區間[1-τ,1]也分為份，但要使得exp(-α(1-x)/ε)等分布在網格節點xi(i=N/2,N/2+1,…,N)上，得到如下網格分布[12]：

顯然，當小參數ε 和網格剖分數N 確定時，網格過渡點τ 由參數α 決定。選取不同的參數α 會使得網格點發生變化，從而影響求解的精度。

在式（2）定義的B-S網格{xi}(i=0,1,…,N)上，記uNi為xi處對應的數值解，可構造問題（1）的迎風差分格式如下：

文獻[12]證明了在B-S網格（2）上利用差分格式（3）求解奇異攝動問題（1）的數值解一致一階收斂于真解，即：

其中，u(xi)和分別為網格點xi處對應的精確解與數值解。

3 優化網格參數的粒子群算法

記uN為當網格剖分為N 個子區間時在B-S網格（2）上求解奇異攝動問題（1）的數值解，它可視為網格參數α 的函數uN(α)。由于很多實際問題中很難得到奇異攝動問題的精確解，因此使用雙重網格的數值解來估計數值解誤差，記：

其中，u2N(α)為將網格剖分為2N 個子區間時的數值解，此時的網格點由網格剖分為N 個子區間時的B-S網格點及相鄰網格的中點組成。

為了使誤差盡可能得小，構造如下的目標函數：

以誤差范數（4）為目標函數尋找最優的網格參數α。

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

其中，ω 為慣性權重，r1,r2是區間（0，1）的隨機數，c1,c2是學習因子。

下面給出采用粒子群算法結合迎風差分格式在Bakhvalov-Shishkin網格上求解奇異攝動對流擴散問題（1）的算法步驟：

步驟1 對參數α 進行種群隨機初始化，由式（2）得到對應的B-S網格。此時迭代次數k=0。

步驟2 利用追趕法等對B-S網格上的差分格式（3）進行求解，得到B-S 網格上N+1 個網格點上的數值解uN(α)。在原B-S 網格上添加相鄰網格點的中點，得到2N+1個網格點，再求解差分格式（3），得到數值解u2N(α)，進而計算出每個粒子對應的目標函數（4）的適應值。

步驟3 對每個粒子α，將其當前適應值與其個體歷史最佳位置(pbest)對應的適應值做比較，如果當前的適應值更優，則將用當前位置更新歷史最佳位置pbest。

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

步驟5 利用式（5）更新粒子α 的速度和位置。由α的新位置，計算對應的Bakhvalov-Shishkin網格。

步驟6 若迭代次數k 小于設定的最大迭代次數Gmax，轉步驟2，并令k=k+1。否則，停止迭代，得到最優網格參數α 及最優誤差范數。

4 數值實驗與結果分析

奇異攝動對流擴散問題（1）中，ε 是擴散項系數，a(x)是對流項系數，當ε 極小時，該問題為穩態的對流占優問題，可以用來描述穩態的熱傳導、粒子擴散、化學反應等現象。下面選擇一個變系數算例和一個常系數算例進行數值模擬。

4.1 變系數奇異攝動問題

考慮奇異攝動問題：

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

4.2 常系數奇異攝動問題

考慮奇異攝動問題：

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

4.3 實驗結果及分析

針對兩個算例，取網格剖分數N =32，64，128，256，小參數ε=10-2,10-3,10-4,10-5,10-6進行實驗。表1 和表2分別列出了以式（4）為目標函數時，算例1和算例2中粒子群算法計算出的最優網格參數及目標函數值?？梢娫谀繕撕瘮底顑灂r，誤差仍保持收斂性。

圖1作出了兩個算例中N=32 ε=10-4時，隨著迭代次數增加目標函數的變化?？梢娊涍^10次左右的迭代，目標函數可接近最優值。針對兩個算例的求解，粒子群算法具有較快的收斂速度。

將參數α 取兩個固定值（分別為搜索范圍的上限和下限）與α 取粒子最優解時的誤差范數目標函數（4）作比較，發現誤差的區別主要體現在x=1 附近的邊界層。問題的解在x=1 附近的邊界層變化劇烈，為了便于觀察邊界層誤差情況，圖2、圖3分別作出了兩個算例中ε=10-4時參數α 取粒子最優解與固定值的邊界層誤差曲線對比，圖（a）和圖（b）分別為網格剖分數N=32 和N=64 的情況。

表1 算例1中粒子群算法計算出的最優網格參數及目標函數值

表2 算例2中粒子群算法計算出的最優網格參數及目標函數值

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

圖2 算例1中ε=10-4 時粒子最優解與固定值的邊界層誤差曲線對比

圖3 算例2中ε=10-4 時粒子最優解與固定值的邊界層誤差曲線對比

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

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

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

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

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

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

5 結束語

本文采用粒子群算法優化選擇網格參數，提出了在Bakhvalov-Shishkin網格上求解一類奇異攝動問題的粒子群算法。實驗結果表明，與選擇固定的網格參數相比，采用粒子群算法計算能得到更好的計算精度，并且數值結果具有收斂性。本文提到的方法也可以推廣到高維問題的計算中。