基于深度學習的偏微分方程求解方法

毛超利

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

0 引言

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

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

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

1 基于深度學習的偏微分方程求解算法介紹

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ò)由若干個神經(jīng)元構(gòu)成。 在整個神經(jīng)網(wǎng)絡(luò)中，如果任意一層神經(jīng)網(wǎng)絡(luò)的任意一個神經(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)的每個神經(jīng)元都包含一個激活函數(shù)（Activation Function）。 激活函數(shù)是簡單的連續(xù)函數(shù)[10]，比如，三角函數(shù)中的雙曲正弦函數(shù)、雙曲正切函數(shù)等。本文選取雙曲正切函數(shù)作為激活函數(shù)，如式（1）所示。

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

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

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

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

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

1.2 本文方法基本框架

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

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

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

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

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

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

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

2 算法驗證

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

2.1 一維對流方程的初值問題

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

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

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

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

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

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

2.2 二維拉普拉斯方程的邊值問題

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

從物理角度， 該方程描述的是一個穩(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é)以第一種提法對算法進行驗證。 給定如下計算域邊界以及邊界條件：

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

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

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

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

2.3 一維擴散方程的初邊值問題

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

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

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

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

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

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

3 結(jié)語

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

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

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