郭曉斌, 袁冬芳, 曹富軍*

(1. 內(nèi)蒙古科技大學 信息工程學院, 內(nèi)蒙古 包頭 014010; 2. 內(nèi)蒙古科技大學 理學院, 內(nèi)蒙古 包頭 014010)

偏微分方程在流體力學、氣候現(xiàn)象、新能源、材料和航天航空領(lǐng)域都有著廣泛的應用[1-2],研究求解偏微分方程的高精度、快速算法具有重要的意義.大多數(shù)偏微分方程都不能找到解析解,因此需要研究數(shù)值近似方法來進行求解.傳統(tǒng)的數(shù)值方法包括有限元法、有限差分法、或有限體積法等[3-5],這些方法需要將復雜計算區(qū)域離散為計算網(wǎng)格,在離散網(wǎng)格上構(gòu)造數(shù)值格式,通過求解數(shù)值格式所形成的代數(shù)方程組獲得方程的近似解.盡管傳統(tǒng)數(shù)值方法在過去的幾十年中已經(jīng)取得了顯著的進步,并且能夠處理相當復雜和高度振蕩的問題,然而數(shù)值解的精度和效率嚴重依賴網(wǎng)格的質(zhì)量.對于復雜計算區(qū)域問題,生成高質(zhì)量的計算網(wǎng)格所花費的時間,在求解時間中占了相當大的比重.此外,因為高維問題中網(wǎng)格點的數(shù)量相對于問題維數(shù)的增加呈指數(shù)級增長,極大增加了問題求解的空間和時間復雜度,不可避免地面臨“維數(shù)災難”等問題.面對以上的挑戰(zhàn), 本文的目標是基于神經(jīng)網(wǎng)絡(luò)的方法求解任意形狀幾何域上定義的橢圓型偏微分方程. 此方法無需生成計算網(wǎng)格, 只需要在計算域內(nèi)部和邊界上進行隨機采樣并作為訓練點集合, 因此可以處理任意形狀的計算區(qū)域,避免了網(wǎng)格生成的困難. 此外,該方法不會受到維數(shù)增加的影響,當維數(shù)增加時,只會增加少量的參數(shù),能夠有效避免“維數(shù)災難”問題,因此具有十分重要的研究意義和價值.

過去幾十年見證了深度學習的革命,人工神經(jīng)網(wǎng)絡(luò)(ANN)的發(fā)展和進化是關(guān)鍵要素之一.基于人工神經(jīng)網(wǎng)絡(luò)的機器學習方法能夠解決圖像處理、語音識別、醫(yī)學診斷等領(lǐng)域一些非常復雜的應用問題. 此外,人工神經(jīng)網(wǎng)絡(luò)在數(shù)學領(lǐng)域也已經(jīng)得到了廣泛的應用,因為它們能夠有效地近似任意函數(shù)[6].雖然人工神經(jīng)網(wǎng)絡(luò)在幾個重要的應用領(lǐng)域取得了令人印象深刻的成果,但關(guān)于它們?yōu)槭裁匆约霸鯓佑行У毓ぷ?仍然存在著許多問題.從傳統(tǒng)意義上來講,使用單一的隱藏層神經(jīng)網(wǎng)絡(luò)可以求解任意的偏微分方程,因為只要一個隱藏層含有足夠多的神經(jīng)元,就可以逼近任何函數(shù),而且神經(jīng)網(wǎng)絡(luò)解的偏導數(shù)可以用解析閉合形式計算[7-8].研究表明,隨著網(wǎng)絡(luò)深度的增加,神經(jīng)網(wǎng)絡(luò)逼近函數(shù)的能力也越強.因此,應用深度神經(jīng)網(wǎng)絡(luò)(DNN)求解偏微分方程得到眾多研究者的廣泛關(guān)注.E等[9]研究了一種求解高維拋物型偏微分方程和倒向隨機微分方程的新算法, 損失函數(shù)由規(guī)定的終端條件和BSDE解之間的誤差給出. E等[10]提出了一種基于深度學習的深度Ritz法,通過變分的方法求解偏微分方程問題.Khoo等[11]提出了一種基于深度神經(jīng)網(wǎng)絡(luò)的求解非齊次偏微分方程的方法,將求解的物理量參數(shù)化,通過神經(jīng)網(wǎng)絡(luò)執(zhí)行時間演化來尋找偏微分方程的解,進而證明使用神經(jīng)網(wǎng)絡(luò)來表征這種量的合理性.Nabian等[12]針對高維隨機偏微分方程提出一種侵入性的、完全不受監(jiān)督的和沒有網(wǎng)絡(luò)的基于深度學習的解決方案.Perdikaris等[13]提出了物理信息神經(jīng)網(wǎng)絡(luò),基于訓練過程實現(xiàn)變分問題的優(yōu)化.Raissi等[14]和Owhadi[15]在機器學習法的基礎(chǔ)上采用高斯過程擬合線性算子,并進一步將該方法推廣到非線性算子的回歸.Berg等[16]基于PINN的方法來求解復雜域上的偏微分方程,并取得了一些良好的結(jié)果.盡管使用PINNs在解決偏微分方程問題上有許多優(yōu)點,但是PINNs仍然存在一些問題.首先,沒有理論基礎(chǔ)來了解神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)的規(guī)模和所需的數(shù)據(jù)量,那么就不能保證算法不會收斂到局部最小值,并且他們的訓練時間要比傳統(tǒng)的數(shù)值方法慢.Dwivedi等[17]開發(fā)了物理信息極限學習機(PIELM),可以應用于和時間相關(guān)的線性偏微分方程問題上,同時證明了PIELM在一系列問題上超過了PINNs的準確性.Anitescu等[18]提出了一種利用人工神經(jīng)網(wǎng)絡(luò)和自適應配置策略求解偏微分方程的方法.

本文研究一種基于深度神經(jīng)網(wǎng)絡(luò)的通用偏微分方程求解方法,結(jié)合物理信息的神經(jīng)網(wǎng)絡(luò)及小樣本學習求解具有復雜計算區(qū)域的問題.該方法不需要生成計算網(wǎng)格或設(shè)計特殊的基函數(shù)或數(shù)值格式來計算橢圓型偏微分方程的解,具有易于實現(xiàn)和無網(wǎng)格的優(yōu)點,有助于求解具有復雜區(qū)域的橢圓型偏微分方程.通過對四類典型的橢圓型偏微分方程進行求解說明了該方法的效率和精度.

1 神經(jīng)網(wǎng)絡(luò)及其結(jié)構(gòu)

首先定義一個單隱層神經(jīng)網(wǎng)絡(luò),給定d維行向量x∈Rd作為模型的輸入,則神經(jīng)網(wǎng)絡(luò)的k維輸出為以下形式

y=σ(xw1+b1)w2+b2

(1)

式中:w1和w2是大小為d×q和q×k的權(quán)重矩陣;b1和b2是大小為和1×k的1×q偏差向量;σ(·)為激活函數(shù).在深度神經(jīng)網(wǎng)絡(luò)中,每一隱層通過權(quán)重矩陣和偏差對輸入變量進行線性處理,然后用激活函數(shù)進行非線性變換. 上一隱層的非線性輸出繼續(xù)作為下一個隱藏層的輸入再次重復權(quán)重矩陣的線性處理及激活函數(shù)的非線性變換.每一隱層中新的權(quán)重矩陣和偏差構(gòu)成了神經(jīng)網(wǎng)絡(luò)中的一個新的隱藏層.通常,神經(jīng)網(wǎng)絡(luò)逼近復雜非線性函數(shù)的能力可以通過增加更多的隱藏層或增加隱藏層的維數(shù)來提升.具有k-1個隱藏層的神經(jīng)網(wǎng)絡(luò)可以表示為

Nθ(z)=Tk°σ°Tk-1°…°T1°σ°T0(z)

(2)

式中:Tk(zk-1)=wkzk-1+bk,wk表示第k-1層的與第k層之間的權(quán)重矩陣,zk-1和bk分別表示第k-1層隱藏層的輸出和偏差向量,Nθ(z)為神經(jīng)網(wǎng)絡(luò)的輸出,z0=z表示輸入?yún)?shù).σ(·)為非線性激活函數(shù).目前,最流行的激活函數(shù)包括Sigmoid、Tanh和ReLU.ReLU激活函數(shù)是使用最廣泛的函數(shù)之一,其形式為f(z)=max(0,z).但是,ReLU的高階導數(shù)為0,這就限制了本文處理高階導數(shù)組成的微分方程的適用性.Tanh或Sigmoid激活函數(shù)可用于二階或高階偏微分方程.Sigmoid激活函數(shù)是非對稱的,將每個神經(jīng)元的輸出限制在區(qū)間[0,1]內(nèi),因此會對神經(jīng)元的輸出引入系統(tǒng)偏差.而Tanh激活函數(shù)是反對稱的,通過允許每個神經(jīng)元的輸出在區(qū)間[-1,1]上取值,克服了Sigmoid激活引起的系統(tǒng)偏差問題.此外,與非對稱激活的深度神經(jīng)網(wǎng)絡(luò)的訓練相比,反對稱激活的深度神經(jīng)網(wǎng)絡(luò)訓練過程收斂的更快.

在給定多個訓練數(shù)據(jù)點的回歸問題中,本文可以使用歐基里德?lián)p失函數(shù)來校準權(quán)重矩陣和偏差,如下所示:

(3)

模型參數(shù)可根據(jù)式(4)進行求解

(4)

使用梯度下降法對模型進行優(yōu)化,具體來說,在第i次迭代中,模型參數(shù)θ={w1,w2,…,b1,b2,…}通過式(5)更新

θ(i+1)=θ(i)-η(i)?θJ(i)(θ(i);x,y)

(5)

式中

η(i)是第i次迭代的步長.

在反向傳播過程中,利用鏈式法則,首先計算最后一層的梯度,最后計算第一層的梯度,當前層的部分梯度在此前的梯度計算中被重復使用.這種信息的反向流動促進了深層神經(jīng)網(wǎng)絡(luò)中每一層梯度的有效計算.算法1是梯度下降迭代算法.該算法主要用于更新神經(jīng)網(wǎng)絡(luò)中的參數(shù).

Algorithm 1：梯度下降迭代算法

1) 給定學習率η;迭代收斂準則eps.

3) while沒有達到停止準則do

采集M個樣本的訓練集{x1,…,xM},對應目標為yi;

計算梯度g=?θJ(θ;x,y);

計算更新:Δθ=-η*g;

應用更新:θ←θ+Δθ;

endwhile

4) 得到神經(jīng)網(wǎng)絡(luò)訓練參數(shù)θ.

2 偏微分方程

根據(jù)以下微分方程求解DNN近似解u(x,θ)

(6)

式中:θ為DNN的可訓練參數(shù);N(·)是一個通用的微分算子,它是由空間導數(shù)以及線性項和非線性項組成;x是有界連續(xù)空間域D?RD上定義的位置向量,D中包含邊界?D.B(·)表示邊界條件,由微分、線性項或非線性項組成.為了計算函數(shù)解,即計算參數(shù)θ,本文通過定義在整個離散域上的非負殘差來計算,

(7)

式中:rB(θ)=rdir(θ)+rneu(θ),rdir(θ)是Dirichlet邊界上的均方誤差;rneu(θ)是Neumman邊界上的均方誤差;Nint和Nbnd分別是計算區(qū)域內(nèi)部點和邊界點的個數(shù);Ndir是Dirichlet邊界點的個數(shù);Nneu是Neumman邊界點的個數(shù),且Nbnd=Ndir+Nneu.

最佳參數(shù)θ*可以根據(jù)式(8)計算,

θ*=argminr(θ)

(8)

因此, 求解微分方程(6)的解便轉(zhuǎn)化為求解式(8)的優(yōu)化問題,其中邊界條件可以視為約束條件.根據(jù)約束條件的處理方法和嚴格程度,將其分為軟約束和硬約束兩種不同的方法.在軟約束中,約束被轉(zhuǎn)化為損失函數(shù)的附加懲罰項[19],這種方法易于實現(xiàn),但不清楚如何調(diào)整損失函數(shù)中不同項的相對重要性,不能保證最終的解決方案會滿足約束條件.在硬約束中,通過構(gòu)造函數(shù),使得任何具有該函數(shù)形式的解都保證滿足邊界條件[7].對于不規(guī)則的邊界,構(gòu)造解函數(shù)的難度很大.本文將采用以上兩種約束方法計算神經(jīng)網(wǎng)絡(luò)近似解.

在軟約束中,令us(x,θ)=uDNN(x,θ)為深度神經(jīng)網(wǎng)絡(luò)的近似解,根據(jù)式(9)求解DNN參數(shù)θ,

(9)

式中:λ1為懲罰項系數(shù),用于調(diào)整損失函數(shù)中每一項的相對重要性[20].

神經(jīng)網(wǎng)絡(luò)通過式(10)更新第i次迭代中的參數(shù)

(10)

(11)

在硬約束中,近似解的函數(shù)形式可以采用以下一般形式[7]:

uh(x,θ)=C(x)+G(x,uDNN(x,θ))

(12)

式中：函數(shù)C(x)滿足邊界條件,函數(shù)G(x,uDNN(x,θ))不受邊界條件的影響.根據(jù)式(13)求解DNN參數(shù)θ,

(13)

為了解決無約束優(yōu)化問題(9)和式(13),本文采用Adam梯度下降算法.在Adam梯度下降算法的每次迭代中,損失函數(shù)的梯度使用輸入空間中的一個點來近似求解,基于該點更新神經(jīng)網(wǎng)絡(luò)的參數(shù).這種迭代更新被證明滿足梯度的無偏估計[22].

同理,在硬約束中,DNN參數(shù)根據(jù)式(14)進行更新,

(14)

式中

(15)

實踐中每次迭代時,損失函數(shù)的梯度是在n個不同的樣本輸入點上計算平均值,而不是僅在一個點上評估.

算法2總結(jié)了神經(jīng)網(wǎng)絡(luò)解法的過程.該算法可以基于預先指定的最大迭代次數(shù)停止,也可以通過損失值的動態(tài)變化作為停止標準.

Algorithm 2：復雜區(qū)域橢圓型偏微分方程的深度學習算法

輸入：區(qū)域內(nèi)部點:xΩ=(x1,…,xi,…,xNint)T,區(qū)域邊界點:xΓ=(x1,…,xj,…,xNbnd)T.

輸出：DNN中的參數(shù)值θ(i).

1) 選擇求解區(qū)域,設(shè)置計算區(qū)域內(nèi)部點和邊界點的個數(shù).

2) 設(shè)置DNN架構(gòu)(層數(shù),每層的維數(shù),激活函數(shù),以及神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)).

3) 初始化DNN參數(shù)θ(0).

4) 設(shè)置Adam優(yōu)化器迭代次數(shù)Niter;設(shè)置學習率η.

5) forifrom 1 toNiterdo

(1) 根據(jù)式(2)計算出神經(jīng)網(wǎng)絡(luò)輸出y;

(3) 通過式(4)和式(5)計算DNN最佳參數(shù)θ(i).

end for

3 數(shù)值算例

本文求解復雜區(qū)域的橢圓型偏微分方程的解,大致步驟為：首先確定計算邊界,并在邊界上選取N個訓練點.同時在計算區(qū)域內(nèi)部任意選取M個訓練點,將兩部分訓練點(M+N)作為訓練集放到tensorflow框架中對神經(jīng)網(wǎng)絡(luò)模型學習,通過訓練得到神經(jīng)網(wǎng)絡(luò)最優(yōu)參數(shù)以及模型的近似解.然后,在計算域內(nèi)選取足夠的點作為測試集,并對訓練后的近似解進行測試,驗證該方法的準確性.

3.1 泊松方程

首先考慮具有混合邊界條件的泊松方程,

其中,x,y∈[0,1].取計算區(qū)域Ω為正方形,在Ω內(nèi)隨機選擇M=500個訓練點,在?Ω上均勻選擇N=200個訓練點,將M和N共同作為訓練集,如圖1所示.且方程的解析解為

圖1 訓練集Fig.1 The traning set

u(x,y)=e-x(x+y3)

損失函數(shù)構(gòu)造如下:

其中,λ1為1 000,λ2為30,選擇學習率η為0.001,經(jīng)過30 000次的迭代訓練,得到方程的訓練解.取M=200 000,N=2 000作為測試集,測試集的解析解和DNN近似解如圖2、圖3所示.

圖2 解析解Fig.2 The analysis solution

圖3 DNN近似解Fig.3 The DNN solution

圖4給出了測試集上的誤差,可以看出神經(jīng)網(wǎng)絡(luò)求解的結(jié)果與解析解的誤差在-1%～5‰之間,效果比較理想.

圖4 測試集誤差Fig.4 The error of the testing sets

針對泊松方程,分別采用Adam,SGD,Adadelta,RMSProp,Adagrad等5種梯度優(yōu)化算法處理損失函數(shù),損失函數(shù)收斂曲線如圖5所示.橫軸代表的是梯度下降算法迭代次數(shù),縱軸表示損失函數(shù)的損失值.隨著迭代次數(shù)的增加,損失值逐漸減小.其中,Adam優(yōu)化算法對模型的優(yōu)化效果要比其余方法要好,主要是因為該算法結(jié)合了Adagrad和RMSProp算法的優(yōu)點,相對于傳統(tǒng)的梯度下降算法來說,Adam優(yōu)化算法具有自適應能力,計算收斂更快.SGD優(yōu)化算法結(jié)果與Adam算法的梯度下降曲線效果相近,但是,當數(shù)據(jù)集比較大的時候,SGD算法的收斂速度和效率優(yōu)于Adam,所以在不同的情況下應該選用合適的優(yōu)化算法.本文計算表明,采用Adam優(yōu)化算法可以達到理想的效果.

圖5 不同優(yōu)化方法的收斂性Fig.5 The convergence of different optimization methods

3.2 二維線形平流方程

平流方程在大氣運動的研究中具有重要意義,是大氣運動中重要的過程之一,運動方程、熱通量方程和水汽方程中都含有平流項.然而,數(shù)值求解方法很難求出此類方程的解析解,因此研究神經(jīng)網(wǎng)絡(luò)解平流方程具有非常重要的現(xiàn)實意義.

其中:x,y∈[0,1],且滿足混合邊界條件,方程的解析解為

u(x,y)=sin(3πx)cos(3πy)

其中,a,b是平流系數(shù),在此問題中,取a=2,b=3,取計算區(qū)域Ω為正方形,在Ω內(nèi)生成M=500個任意的訓練點,在?Ω上生成N=200個訓練點,如圖6所示.

圖6 訓練集Fig.6 The traning set

圖7、圖8分別為測試集在解析解和DNN解上的函數(shù)結(jié)果,圖9給出了二維線性平流方程在測試集上的誤差,誤差結(jié)果在-4‰～8‰之間.

圖7 解析解Fig.7 The analysis solution

圖8 DNN近似解Fig.8 The DNN solution

圖9 測試集誤差Fig.9 The error of the testing set

損失函數(shù)構(gòu)造如下:

其中:λ1=200,λ2=50,λ3=50,Nneu1是y=1時的Neumman邊界點的個數(shù),Nneu2是y=0時的Neumman邊界點的個數(shù).取M=200 000,N=2 000作為測試集.

圖10比較了隨著迭代次數(shù)的增加,不同優(yōu)化算法的收斂曲線.

圖10 不同優(yōu)化方法的收斂性Fig.10 The convergence of different optimization methods

3.3 二維對流擴散方程

本算例考慮二維對流擴散方程

其中:計算區(qū)域為花瓣形,花瓣的半徑[23]可以表示為Φ(x,y)=r-(0.4+0.3sin(5θ)),其中

將Ω區(qū)間劃分為80×80個網(wǎng)格點,判斷網(wǎng)格點是否落在花瓣半徑以內(nèi),當-0.1<Φ(xi,yi)<0時,此時的(xi,yi)i=1,2,…,M是花瓣的邊界,當Φ(xj,yj)<-0.1,(xj,yj)j=1,2,…,N為花瓣的內(nèi)部點,如圖11所示.

圖11 訓練集Fig.11 The traning set

方程的解析解為

損失函數(shù)構(gòu)造如下:

其中:λ1為1 000選取花瓣內(nèi)的網(wǎng)格點M=14 400,N=6 300作為測試集.圖12和圖13給出了測試集在解析解和DNN近似解函數(shù)圖像.圖14給出了測試集的誤差,誤差范圍在-6‰～4‰之間,表現(xiàn)出較高的精度.圖15比較了隨著迭代次數(shù)的增加,不同優(yōu)化方法的損失函數(shù)收斂曲線,可以看出Adam方法的收斂效率優(yōu)于其它四種方法.

圖12 解析解Fig.12 The analysis solution

圖13 DNN近似解Fig.13 The DNN solution

圖14 測試集誤差Fig.14 The error of the testing set

圖15 不同優(yōu)化方法的收斂性Fig.15 The convergence of different optimization methods

3.4 二維勢能方程——由梯度勢能計算力場

在分子動力學中,力場可以由勢能的負梯度求得

方程的解析解為

u(x,y)=1+y2-2x2+x6

損失函數(shù)構(gòu)造如下：

其中：λ1為100,F(x,y)是力場,u(x,y)是勢能.計算區(qū)域Ω為中國大陸,在Ω內(nèi)任意選擇M=500個訓練點,在?Ω上選擇N=455個訓練點,如圖16所示.采用具有2個隱藏層、30個神經(jīng)元的神經(jīng)網(wǎng)絡(luò)計算訓練解.取M=100 000,N=2 276作為測試集代入ANN訓練解中.圖17為測試集上的精確解.圖18為測試集上的近似解.圖19給出了測試集上的誤差,誤差范圍在-1%～7‰之間.圖20比較了隨著迭代次數(shù)的增加,不同優(yōu)化方法的損失收斂曲線.由于不同偏微分方程的復雜度是不同的,復雜度高的偏微分方程在訓練過程中不容易訓練,損失函數(shù)收斂過程比較慢,且不易收斂,在訓練過程中需要增加訓練次數(shù),并采用適當?shù)膬?yōu)化方法來處理.

圖16 訓練集Fig.16 The traning set

圖17 解析解Fig.17 The analysis solution

圖18 DNN近似解Fig.18 The DNN solution

圖19 測試集誤差Fig.19 The error of the testing set

圖20 不同優(yōu)化方法的收斂性Fig.20 The convergence of different optimization methods

4 結(jié)論

本文研究一種利用神經(jīng)網(wǎng)絡(luò)求解復雜區(qū)域橢圓型偏微分方程的方法.該方法使用完全無監(jiān)督的類梯度下降算法進行迭代訓練,在求解過程中無需調(diào)用其他的偏微分方程數(shù)值解算器.此外,它在整個空間域上是無網(wǎng)格的,因此,有利于處理不規(guī)則的計算區(qū)域.通過反向傳播算法的推導,將二階偏微分方程擴展到任意階偏微分方程組.數(shù)值算例結(jié)果表明,在本文提供的神經(jīng)網(wǎng)絡(luò)計算框架下,任意橢圓型偏微分方程在任何區(qū)域內(nèi)都可以通過訓練得到它們的近似解,且具有較高的精度;同時,神經(jīng)網(wǎng)絡(luò)方法求解復雜區(qū)域橢圓型偏微分方程的收斂性和計算效率還需要進一步提升和改進.因此,還需要考察不同的優(yōu)化方法,進一步提高神經(jīng)網(wǎng)絡(luò)方法的性能.

致謝：本文得到內(nèi)蒙古科技大學創(chuàng)新基金(2019YQL02)項目的資助,在此表示感謝.