毛超利

（中國(guó)電子科技南湖研究院，浙江 嘉興 314002）

0 引言

伴隨著計(jì)算機(jī)技術(shù)的高速發(fā)展，科學(xué)計(jì)算已經(jīng)成為可與理論和實(shí)驗(yàn)科學(xué)方法相并列的第三種科學(xué)方法[1][2]。科學(xué)計(jì)算中的一個(gè)重要內(nèi)容是數(shù)值求解各類偏微分方程（組）[3]。 數(shù)值求解偏微分方程（組）具有重要意義。 比如，核爆炸物理過(guò)程和流體宏觀流動(dòng)均可用偏微分方程組來(lái)描述，借助理論手段解析地求解這些偏微分方程，在當(dāng)今數(shù)學(xué)發(fā)展水平下依然不可能[4]，利用科學(xué)計(jì)算方法求取它們的數(shù)值逼近解，可以補(bǔ)充試驗(yàn)無(wú)法得到的細(xì)節(jié)，從而利于減少試驗(yàn)次數(shù)、節(jié)約研發(fā)成本和縮短研發(fā)周期。 有資料表明[1]，從開始研制至原子彈爆炸成功，我國(guó)進(jìn)行了338 次核試驗(yàn)，而美國(guó)、蘇聯(lián)分別進(jìn)行了936次和716 次試驗(yàn)，其中一個(gè)重要原因是我國(guó)科學(xué)家充分采用了基于科學(xué)計(jì)算方法求解偏微分方程的數(shù)值模擬手段。

數(shù)值求解偏微分方程的算法包括有限差分法、有限體積法、有限元法和譜方法。 有限差分法是用差商近似代替微分[2]，從而將連續(xù)的偏微分方程轉(zhuǎn)化為代數(shù)方程進(jìn)行求解。有限體積法是將空間離散為結(jié)構(gòu)化網(wǎng)格[2]，在每個(gè)網(wǎng)格內(nèi)施加守恒定律，從而將全局的偏微分方程轉(zhuǎn)化為局部代數(shù)方程進(jìn)行求解。有限元法是把原來(lái)待求解的幾何空間劃分為有限元，并假設(shè)每個(gè)有限元上的原函數(shù)關(guān)系可以近似地由一組基函數(shù)進(jìn)行線性組合，將該線性組合帶入原偏微分方程，會(huì)有一個(gè)誤差余量，通過(guò)引入加權(quán)函數(shù)將誤差余量最小化，就可以得到一個(gè)線性方程組，對(duì)其求解即可得到原偏微分方程的數(shù)值解[5]。譜方法與有限元法類似，不同之處在于，譜方法是在全局時(shí)空上構(gòu)造原函數(shù)基于基函數(shù)的線性組合[6]。以上四種方法均已得到充足的發(fā)展。 但是，它們各自有著一些無(wú)法克服的缺點(diǎn)。 比如，有限差分法會(huì)引入人工粘性。

近些年來(lái)，隨著深度學(xué)習(xí)在眾多實(shí)際應(yīng)用中取得成功，越來(lái)越多的研究人員開始嘗試?yán)盟鼇?lái)解決各自領(lǐng)域的傳統(tǒng)難題[7-9]。 受此啟發(fā)，針對(duì)偏微分方程數(shù)值求解問(wèn)題，本文嘗試了一種新思路——基于深度學(xué)習(xí)的方法。該方法與有限元法和譜方法類似，都是對(duì)偏微分方程真實(shí)解對(duì)應(yīng)非線性關(guān)系的一種數(shù)值逼近。 不同之處在于，有限元方法和譜方法是一次輸入對(duì)應(yīng)一次輸出的單一映射，而在本文方法中，一層神經(jīng)網(wǎng)絡(luò)的輸出又是下一層神經(jīng)網(wǎng)絡(luò)的輸入，因此具有更強(qiáng)的非線性擬合能力。

1 基于深度學(xué)習(xí)的偏微分方程求解算法介紹

1.1 深度神經(jīng)網(wǎng)絡(luò)

前饋深度神經(jīng)網(wǎng)絡(luò)是最基礎(chǔ)的神經(jīng)網(wǎng)絡(luò)，由輸入層、隱藏層和輸出層構(gòu)成，信息由輸入層向輸出層單向傳播。輸入層和輸出層一般具有單層神經(jīng)網(wǎng)絡(luò)，隱藏層可具有多層神經(jīng)網(wǎng)絡(luò)，每一層神經(jīng)網(wǎng)絡(luò)由若干個(gè)神經(jīng)元構(gòu)成。 在整個(gè)神經(jīng)網(wǎng)絡(luò)中，如果任意一層神經(jīng)網(wǎng)絡(luò)的任意一個(gè)神經(jīng)元和下一層神經(jīng)網(wǎng)絡(luò)的所有神經(jīng)元都是連接的，則稱這種神經(jīng)網(wǎng)絡(luò)是全連接的[10]。 本文使用的就是這種前饋全連接神經(jīng)網(wǎng)絡(luò)，如圖1 所示。

圖1 神經(jīng)網(wǎng)絡(luò)示意圖Figure 1 Schematic diagram of neural network

隱藏層內(nèi)的每個(gè)神經(jīng)元都包含一個(gè)激活函數(shù)（Activation Function）。 激活函數(shù)是簡(jiǎn)單的連續(xù)函數(shù)[10]，比如，三角函數(shù)中的雙曲正弦函數(shù)、雙曲正切函數(shù)等。本文選取雙曲正切函數(shù)作為激活函數(shù)，如式（1）所示。

通常，隱藏層內(nèi)的所有神經(jīng)元包含同一形式的激活函數(shù)。將與隱藏層某一神經(jīng)元連接的其他神經(jīng)元傳輸進(jìn)來(lái)的信息進(jìn)行加權(quán)（weighted）平均，再加上一個(gè)偏置（bias），作為該神經(jīng)元的輸入，在激活函數(shù)的作用下，該神經(jīng)元輸出信息，并傳遞給下一層神經(jīng)元。

用矩陣的形式表達(dá)，隱藏層的第一層接收到的信息如式（2）所示。

式（3）分別代表輸入層作用于隱藏層第一層的權(quán)重系數(shù)矩陣和輸入信息向量。在激活函數(shù)的作用下，隱藏層第一層神經(jīng)元的輸出信息如式（4）所示。

其中：σ 代表激活函數(shù)，B1代表該層神經(jīng)元的偏置向量。對(duì)于圖1 所示包含兩個(gè)隱藏層的神經(jīng)網(wǎng)絡(luò)，信息經(jīng)過(guò)過(guò)濾，到達(dá)輸出層時(shí)，如式（5）所示。

其中，W3∈R4×1是隱藏層第二層向輸出層傳遞時(shí)的權(quán)重系數(shù)矩陣。 可以看到，最終的輸出信息是關(guān)于輸入信息和一系列權(quán)重系數(shù)以及偏置的復(fù)雜非線性關(guān)系。 選取合適的權(quán)重系數(shù)和偏置，則由神經(jīng)網(wǎng)絡(luò)表示的非線性關(guān)系是對(duì)原偏微分方程解的一個(gè)近似。

1.2 本文方法基本框架

為了討論的方便，不失一般性地，本文考慮如下形式的二階偏微分方程，如式（6）所示。

其中，C 代表任意常數(shù)項(xiàng)。 給定初始條件，或邊界條件，或混合條件（初始條件加上邊界條件），偏微分方程的解是唯一確定的。 其次，假設(shè)使用圖1 所示的神經(jīng)網(wǎng)絡(luò)對(duì)上述偏微分方程的解進(jìn)行逼近，則：

其中，右邊項(xiàng)中的參數(shù)符號(hào)意義與第1.1 節(jié)保持一致，a=[t，x，y，z]T是輸入信息向量。 一方面，由于uNN不是偏微分方程的精確解，故將其代入原偏微分方程左端時(shí)，結(jié)果不等于0，而是有一個(gè)余量偏差 Rf，如式（8）所示。

這里需要指出的是，在深度神經(jīng)網(wǎng)絡(luò)中，輸出信息關(guān)于輸入信息的偏微分可以使用自動(dòng)微分（Auto-difference）計(jì)算，自動(dòng)微分本質(zhì)上是求導(dǎo)鏈?zhǔn)椒▌t，具有很高的效率。另一方面，對(duì)應(yīng)初始時(shí)刻的輸出信息與初始條件的余量偏差Ri為：

在求解域上均勻分布的一系列時(shí)空坐標(biāo)下，取值最小。 其中，Nx，Ny，Nz和 Nt分別是計(jì)算域在 x，y，z 和 t 方向上的離散點(diǎn)個(gè)數(shù)，NxΓ，NyΓ，NzΓ分別代表計(jì)算域在x，y，z 邊界上的離散點(diǎn)個(gè)數(shù)。 選取偏差余量的平方和作為準(zhǔn)則的依據(jù)與最小二乘法的依據(jù)類似，具體細(xì)節(jié)可以參考相關(guān)文獻(xiàn)[5]。 式（11）即是基于深度學(xué)習(xí)偏微分方程求解算法的最小化目標(biāo)函數(shù)， 它可以看作是關(guān)于一系列權(quán)重系數(shù)W 和偏置B 的函數(shù)，即：

因此，偏微分方程求解問(wèn)題最終轉(zhuǎn)化為一個(gè)無(wú)約束最優(yōu)化問(wèn)題。求解無(wú)約束最優(yōu)化問(wèn)題常用的數(shù)值算法有梯度下降法[11]、擬牛頓算法BFGS[12]等。 擬牛頓算法BFGS 本質(zhì)上是從梯度下降法開始迭代，不斷逼近牛頓法求曲線駐點(diǎn)的算法，其效率比梯度下降法高。 經(jīng)過(guò)研究人員的不斷改進(jìn)，L-BFGS 算法[12]進(jìn)一步克服了原始的BFGS 算法內(nèi)存需求大的缺點(diǎn)。 因此，本文采用L-BFGS 算法求解由偏微分方程求解轉(zhuǎn)化而來(lái)的無(wú)約束最優(yōu)化問(wèn)題。 綜上所述，本文基于深度學(xué)習(xí)的偏微分方程求解算法基本框架如圖2 所示。

圖2 算法流程圖Figure 2 Algorithm flow chart

2 算法驗(yàn)證

從數(shù)學(xué)角度， 偏微分方程可以分為橢圓型方程、雙曲型方程和拋物型方程[3]。本文選取對(duì)應(yīng)這種分類的三個(gè)算例對(duì)算法進(jìn)行驗(yàn)證。

2.1 一維對(duì)流方程的初值問(wèn)題

一維對(duì)流方程是典型的一階雙曲型偏微分方程，具有如下形式：

其中，a 是非零常數(shù)。 我們知道， 它是雙曲型的，給定初始條件，它的解是唯一的。

不失一般性地，這里令a=-1，并給定如下初始條件：

初始信號(hào)是三角形的，寬度為0.2，且在x=0 處一階導(dǎo)數(shù)不連續(xù)。 為了檢驗(yàn)本文算法準(zhǔn)確性，分別采用本文算法和有限差分迎風(fēng)格式算法求解上述雙曲型偏微分方程的初值問(wèn)題。不同時(shí)刻的結(jié)果如圖3 所示。

圖3 采用迎風(fēng)格式算法和本文深度學(xué)習(xí)算法計(jì)算結(jié)果對(duì)比Figure 3 Comparison of calculation results using upwind style and deep learning algorithm

理想情況下，求解的結(jié)果應(yīng)該是在不同時(shí)刻分布在不同空間位置的三角形信號(hào)形狀。從圖3 可以看到，由于引入了“人工粘性”[2]，迎風(fēng)格式“抹平”了信號(hào)尖角，而且隨著時(shí)間的演化，信號(hào)強(qiáng)度也被削弱。 相比之下，本文基于深度學(xué)習(xí)算法的結(jié)果較好地保持了初始信號(hào)的形狀和強(qiáng)度。

2.2 二維拉普拉斯方程的邊值問(wèn)題

二維拉普拉斯方程是典型的橢圓型偏微分方程，不失一般性地，這里假設(shè)它具有如下形式：

從物理角度， 該方程描述的是一個(gè)穩(wěn)態(tài)現(xiàn)象，比如在邊界溫度給定的條件下，一塊金屬平板內(nèi)部溫度的分布。橢圓型偏微分方程的邊界條件有以下三種提法：（1）固定邊界條件，即在給定邊界上給定u 的值 U1（x，y）；（2）在給定邊界上給定 u 的法向?qū)?shù)值，=U2（x，y）；（3）在給定邊界上給定混合邊界條件，+u=U3（x，y）。 其中，第一種提法最為普遍， 本小節(jié)以第一種提法對(duì)算法進(jìn)行驗(yàn)證。 給定如下計(jì)算域邊界以及邊界條件：

分別采用本文算法和五點(diǎn)差分格式求解。整體分布云圖定性對(duì)比和不同位置處的結(jié)果定量對(duì)比如圖4 所示。 可以看到，本文算法與五點(diǎn)差分格式具有相當(dāng)?shù)木取?/p>

圖4 二維拉普拉斯方程數(shù)值解結(jié)果Figure 4 Numerical solution results of the two-dimensional Laplace equation

更進(jìn)一步地，本小節(jié)探討深度神經(jīng)網(wǎng)絡(luò)隱藏層層數(shù)和寬度（每層的神經(jīng)元數(shù)目）對(duì)計(jì)算結(jié)果的影響。 相對(duì)誤差百分比的計(jì)算選取五點(diǎn)差分格式結(jié)果作為精確解。 結(jié)果如表1 所示。 可以看出，對(duì)于二維拉普拉斯方程，隨著隱藏層寬度的增加，神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)精度大體上先上升后下降， 而隨著隱藏層層數(shù)的增加， 預(yù)測(cè)精度大體上反而下降。 這說(shuō)明， 使用神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)二維拉普拉斯方程邊值問(wèn)題數(shù)值解時(shí)， 采用更多的神經(jīng)元并不一定會(huì)帶來(lái)預(yù)測(cè)精度的提升。

表1 二維拉普拉斯方程神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)數(shù)值解的相對(duì)誤差百分比Table 1 The relative error percentage of the twodimensional Laplace equation neural network predicting the numerical solution

2.3 一維擴(kuò)散方程的初邊值問(wèn)題

擴(kuò)散方程是應(yīng)用中常見(jiàn)的拋物型偏微分方程，它有著明顯的物理背景，例如不同濃度物質(zhì)之間的擴(kuò)散過(guò)程、熱傳遞過(guò)程和波傳播過(guò)程等。 如式（17）所示的偏微分方程：

式（17）是典型形式的擴(kuò)散方程。 其中，a 是非零常數(shù)。 給定初始和邊界條件，上述偏微分方程有唯一解。

選取a=-1，并給定如下初始以及邊界條件：

分別采用本文算法和Crank-Nicolson 格式求解。神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)的整體分布云圖和不同時(shí)刻的結(jié)果定量對(duì)比如圖5 所示。 其中圖5a）為神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)值云圖，白線標(biāo)記定量對(duì)比的時(shí)刻。可以看到，本文算法與Crank-Nicolson 格式具有相當(dāng)?shù)木取?/p>

圖5 一維擴(kuò)散方程數(shù)值解結(jié)果Figure 5 Numerical solution results of one-dimensional diffusion equation

更進(jìn)一步地，本小節(jié)探討深度神經(jīng)網(wǎng)絡(luò)隱藏層層數(shù)和寬度（每層的神經(jīng)元數(shù)目）對(duì)計(jì)算結(jié)果的影響。相對(duì)誤差百分比的計(jì)算選取Crank-Nicolson 格式結(jié)果作為精確解。結(jié)果如表2 所示?？梢钥闯觯愃朴诙S拉普拉斯方程結(jié)果， 對(duì)于一維擴(kuò)散方程，隨著隱藏層寬度的增加，神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)精度大體上先上升后下降，而隨著隱藏層層數(shù)的增加，預(yù)測(cè)精度大體上反而下降。 這說(shuō)明，使用神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)一維擴(kuò)散方程初邊值問(wèn)題數(shù)值解時(shí)，采用更多的神經(jīng)元并不一定會(huì)帶來(lái)預(yù)測(cè)精度的提升。

3 結(jié)語(yǔ)

本文提出了一種基于深度學(xué)習(xí)的偏微分方程求解方法。 針對(duì)雙曲型、橢圓型和拋物型偏微分方程三種典型問(wèn)題，分別使用有限差分格式和本文方法進(jìn)行求解。 結(jié)果對(duì)比表明，本文算法具有比較好的求解精度，而且它能夠克服迎風(fēng)格式引入人工粘性的缺點(diǎn)。 此外，本文算法不需要針對(duì)不同類型的偏微分方程作特殊處理，即本文算法框架具有更廣泛的適用性。

表2 一維擴(kuò)散方程神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)數(shù)值解的相對(duì)誤差百分比Table 2 The relative error percentage of the onedimensional diffusion equation neural network predicting the numerical solution

本文方法的以上優(yōu)點(diǎn)并不意味著目前它可以在工程應(yīng)用中取代業(yè)已發(fā)展成熟的傳統(tǒng)方法。原因是，基于深度神經(jīng)網(wǎng)絡(luò)求解偏微分方程時(shí)，誤差規(guī)律還不清楚。 此外，使用深度學(xué)習(xí)方法求解偏微分方程要求問(wèn)題的邊界條件或初始條件有比較明確的定義，否則可能會(huì)出現(xiàn)邊界捕捉不到或者瞬態(tài)求解不準(zhǔn)確的問(wèn)題。鑒于以上兩個(gè)方面，未來(lái)，我們需要對(duì)訓(xùn)練過(guò)程中深度神經(jīng)網(wǎng)絡(luò)是如何收斂到偏微分方程解的數(shù)學(xué)本質(zhì)進(jìn)行深入研究。