趙青宇， 鄭素佩， 李 霄

(長安大學 理學院，西安 710064)

1 引 言

雙曲守恒律方程是描述空氣動力學、物理學和海洋學等諸多領域問題的一類數(shù)學模型?；诓徽摮跏紬l件多么光滑，該類方程的解都可能會出現(xiàn)間斷這一事實，其數(shù)值求解算法一直廣受關注。針對非線性雙曲型守恒律方程的間斷現(xiàn)象，經(jīng)典解理論不再適用，于是引入了弱解的概念，并提出通過熵穩(wěn)定條件來保證解的唯一性。為滿足這一條件，應添加適當?shù)臄?shù)值粘性項。而對于這一問題，Tadmor[1]定義了一類具有二階精度的熵守恒格式。為進一步研究，根據(jù)Tadmor[1]給出的比較原則，Ismail等[2]提出了熵穩(wěn)定格式。Fjordholm等[3]基于滿足符號性質的ENO[4]重構了高階熵穩(wěn)定格式。鄭素佩等[5]在熵穩(wěn)定格式中利用不同高階重構得到高精度熵穩(wěn)定格式。呂夢迪等[6]采用五階CWENO[7]重構得到一種新的熵穩(wěn)定算法。胡立軍等[8]采用通量分裂方法對此類方程求解。除此之外，網(wǎng)格的分布對數(shù)值結果的好壞有著較大的影響，Harten等[9]對移動網(wǎng)格法的發(fā)展做出了較大貢獻，這一方法優(yōu)勢在于對網(wǎng)格的有效利用，在物理解變化較大的區(qū)域網(wǎng)格自動加密，反之網(wǎng)格自動稀疏。到目前為止，該方法的發(fā)展已取得較多成果。楊繼明[10]對求解奇異攝動問題的移動網(wǎng)格算法進行了研究。根據(jù)Tang等[11]提出的自適應移動網(wǎng)格法,文獻[12，13]基于移動網(wǎng)格的熵穩(wěn)定格式對雙曲型守恒律方程進行了求解。

近年來，神經(jīng)網(wǎng)絡得到了快速發(fā)展，Lagaris等[14]提出了一種用人工神經(jīng)網(wǎng)絡求解初邊值問題的方法；Yohai等[15]提出一種基于已有方程的數(shù)值解對偏微分方程優(yōu)化逼近的方法；張琪[16]提出了一種訓練基于物理知識的神經(jīng)網(wǎng)絡求解偏微分方程的新方法。

傳統(tǒng)數(shù)值算法在間斷附近，結果常會存在偽振蕩或過度抹平現(xiàn)象。本文基于機器學習框架下的 BP神經(jīng)網(wǎng)絡訓練現(xiàn)有數(shù)值方法中的函數(shù)值，給出一種求解雙曲守恒律方程數(shù)值解的新算法，用于提高數(shù)值結果的分辨率和避免數(shù)值偽振蕩的產(chǎn)生，數(shù)值結果表明新算法具有魯棒性強和分辨率高的特點。

2 BP神經(jīng)網(wǎng)絡

2.1 算法描述

BP神經(jīng)網(wǎng)絡是一種根據(jù)誤差反向傳播訓練的多層前饋神經(jīng)網(wǎng)絡，其基本思想[17]為，學習過程由正向傳播與反向傳播組成，根據(jù)隱含層每個神經(jīng)元的輸入與輸出，得到輸出層的輸入與輸出，再將網(wǎng)絡輸出與期望輸出進行對比，若正向傳播的輸出結果誤差達到理想誤差，則停止計算，否則利用梯度下降法優(yōu)化誤差，計算出各層的權重與閾值的修正項后繼續(xù)反向傳播，重復迭代。

如圖1所示，一個三層的BP神經(jīng)網(wǎng)絡，用i,j和k分別表示各層數(shù)據(jù)節(jié)點。學習過程如下。

圖1 BP神經(jīng)網(wǎng)絡結構

(1) 正向傳播

設隱含層的輸入數(shù)據(jù)和輸出數(shù)據(jù)分別為Hinput和Houtput，輸入層和隱含層權值為Wi j，隱含層與輸出層的權值為Wj k，隱含層和輸出層的閾值分別為b1和b2,激活函數(shù)為f，則

Hinput j=∑Wi jxi-b1j,Houtput j=f(Hinput j)

設輸出層的輸入數(shù)據(jù)和輸出數(shù)據(jù)分別為E和Y,則

Ek=∑Wj kHoutput j-b2k,Yk=f(Ek)

網(wǎng)絡誤差D表示為

(2) 反向傳播

對于輸出層和隱含層的節(jié)點，修正項分別為

δk=(Ak-Yk)Yk(1-Yk)

δj=Houtput j(1-Houtput j)∑δkWjk

反向修正后輸入層與隱含層的權值和閾值為

對上述過程進行循環(huán)迭代，直至達到理想誤差。

2.2 數(shù)據(jù)集的構造

將熵穩(wěn)定格式和基于移動網(wǎng)格的熵穩(wěn)定格式作為已知數(shù)值格式(具體在第3節(jié)給出)，采用已知格式的低分辨率網(wǎng)格點數(shù)據(jù)作為網(wǎng)絡輸入，其中一種高分辨率網(wǎng)格點數(shù)據(jù)作為期望輸出，構造訓練集時，用時間方向的多層輸入數(shù)據(jù)表示。構造預測集時，用任一給定時間層的網(wǎng)格點數(shù)據(jù)表示(具體構造在算例中說明)。

2.3 網(wǎng)絡模型

BP神經(jīng)網(wǎng)絡求解雙曲守恒律方程時，使用均方誤差定義損失函數(shù)，網(wǎng)絡訓練時均采用trainlm函數(shù)，網(wǎng)絡模型包含兩個隱層對輸入數(shù)據(jù)進行非線性處理，在求解中設置隱層神經(jīng)元個數(shù)分別為23和18，激活函數(shù)為雙曲正切函數(shù)，其表達式為

tanhx=(ex-e- x)/(ex+e- x)

網(wǎng)絡訓練與預測流程如圖2所示。

圖2 網(wǎng)絡訓練與預測流程

3 數(shù)值算法

3.1 熵穩(wěn)定格式

一維雙曲守恒律方程

Ut+F(U)x=0

(1)

式中U=U(x,t)∈Rm為守恒變量,F(U)為通量,式(1)的半離散格式為

(2)

式中Δxi=xi + 1/2-xi - 1/2,Ui(t)為U(x,t)在網(wǎng)格單元[xi - 1/2,xi + 1/2]的單元平均,Fi + 1/2為與F(U)相容的數(shù)值通量。方程(1)的熵穩(wěn)定格式參見文獻[1]。

3.2 移動網(wǎng)格法

定義物理區(qū)域Ωp和邏輯區(qū)域Ωc存在一一對應的映射,即x=x(ξ),ξ∈Ωc,x∈Ωp，首先，對式(2)中心差分格式離散化，并采用Gauss-Seidel迭代求解方程實現(xiàn)網(wǎng)格的更新

其次，采用守恒型插值公式實現(xiàn)解的更新

最后，為減少數(shù)值解誤差，選取如下低通濾波器使網(wǎng)格函數(shù)更光滑，參見文獻[11]。

ω←0.5ωi+0.25(ωi - 1-ωi + 1)

4 數(shù)值算例

本節(jié)先說明數(shù)據(jù)的選擇，再通過不同數(shù)值算例結果來驗證新方法的性能。

假設方程(1)取N個網(wǎng)格點計算時間T，每個時間層均選用熵穩(wěn)定格式和移動網(wǎng)格法計算節(jié)點處函數(shù)值，將T的前m個時間層所得的各個節(jié)點處計算函數(shù)值作為網(wǎng)絡輸入，將對應時間層的高分辨率數(shù)值格式的計算函數(shù)值作為網(wǎng)絡輸出，其中這兩組數(shù)據(jù)集作為訓練集。同樣將該時間層T采用熵穩(wěn)定格式和移動網(wǎng)格法的計算函數(shù)值作為預測集。

算例1考慮一維線性對流方程滿足間斷初始條件

取100個網(wǎng)格點計算時間t=0.5,CFL條件數(shù)為0.25。該算例的準確解含有激波，圖3(a,b)分別給出了熵穩(wěn)定格式和移動網(wǎng)格法的結果比較，圖3(c)是由50個時間層的數(shù)據(jù)訓練結果與當前時間層數(shù)據(jù)仿真預測的結果，可以看出，新方法在間斷位置計算結果明顯比已有算法好，神經(jīng)網(wǎng)絡能很好地學習到方程的結構并捕捉到激波，與參考解非常接近。圖3(d)給出了網(wǎng)絡收斂性能，若目標精度為1×10-5，則神經(jīng)網(wǎng)絡在迭代393次后誤差精度達到5.2257×10-6，性能良好。

圖3 算例1的數(shù)值結果與損失函數(shù)變化

算例2考慮一維無粘Burgers方程滿足間斷初始條件

ut+(u2/2)x=0

取100個網(wǎng)格點計算時間t=0.3,CFL條件數(shù)為0.2。該算例的準確解是由稀疏波和激波兩部分組成，圖4(a,b)分別給出了熵穩(wěn)定格式和移動網(wǎng)格法的結果比較，圖4(c)是由40個時間層的數(shù)據(jù)訓練結果與當前時間層數(shù)據(jù)仿真所預測的結果，可以看出，新方法在激波和稀疏波間斷位置的結果與參考解誤差很小，對已有格式的抹平現(xiàn)象明顯改善。圖4(d)給出了網(wǎng)絡收斂性能，若目標精度為1×10-5，則神經(jīng)網(wǎng)絡在迭代232次后誤差精度達到8.6412×10-6，性能良好。

圖4 算例2的數(shù)值結果與損失函數(shù)變化

考慮一維無粘Burgers方程壓縮波滿足間斷初始條件

取100個網(wǎng)格點計算時間t=0.96,CFL條件數(shù)為0.2，采取周期邊界條件，該問題的準確解含有間斷，產(chǎn)生了激波，圖5(a,b)分別給出了熵穩(wěn)定格式和移動熵穩(wěn)定格式結果比較，圖5(c)是由40個時間層的數(shù)據(jù)訓練結果與當前時間層數(shù)據(jù)仿真預測的結果，可以看出，該方法的計算結果既保持已有算法在光滑區(qū)域的良好數(shù)值解，又可以根據(jù)多個時間層數(shù)據(jù)學習到方程的結構，在間斷位置比熵穩(wěn)定格式和移動熵穩(wěn)定格式更接近于參考解，結果更準確。圖5(d)給出了網(wǎng)絡收斂性能，若目標精度為1×10-5，則神經(jīng)網(wǎng)絡在迭代211次后誤差精度達到9.857×10-6，性能良好。

算例4一維帶源項淺水波方程在間斷底部上的潰壩問題

計算區(qū)間[0,150]，初始條件和底部地勢函數(shù)分別為

u(x,0)=0

取200個網(wǎng)格點計算時間t=4.5,CFL條件數(shù)為0.1,g=0.98，采用周期性邊界條件。該算例的準確解是由稀疏波和激波兩部分組成，可以看出，圖6(a,b)分別給出了熵穩(wěn)定格式和移動網(wǎng)格法的結果比較，圖6(c)是由40個時間層的數(shù)據(jù)訓練結果與當前時間層數(shù)據(jù)仿真預測的結果，在間斷位置預測解可以自動識別，結果比熵穩(wěn)定格式和移動熵穩(wěn)定格式更貼近參考解。圖6(d)給出了網(wǎng)絡收斂性能，若目標精度為 1×10-5，則神經(jīng)網(wǎng)絡在迭代324次后誤差精度達到5.3882×10-5，性能良好。

圖5 算例3的數(shù)值結果與損失函數(shù)變化

圖6 算例4的數(shù)值結果與損失函數(shù)變化

算例5一維帶源項淺水波方程捕捉穩(wěn)態(tài)水中出現(xiàn)的微小擾動問題

計算區(qū)間[0,2]，初始條件和連續(xù)形函數(shù)分別為

hu(x,0)=0

取200個網(wǎng)格點計算時間t=0.3,CFL條件數(shù)為0.45,其中g=0.98，ξ=0.001，采用周期性邊界條件。該算例是一個連續(xù)問題，圖7(a,b)分別給出了熵穩(wěn)定格式和移動網(wǎng)格法的結果比較，圖7(c)是由40個時間層的數(shù)據(jù)訓練結果與當前時間層數(shù)據(jù)仿真預測的結果，可以看出采用新方法求解連續(xù)問題時，在斜率變化較大和斜率為零的位置新方法計算結果更好，更接近參考解。圖7(d)給出了網(wǎng)絡收斂性能，若目標精度為1×10-5，則神經(jīng)網(wǎng)絡在迭代306次后誤差精度達到9.9931×10-6，性能良好。

圖7 算例5的數(shù)值結果與損失函數(shù)變化

5 結 論

本文采用熵穩(wěn)定格式和移動熵穩(wěn)定格式兩種經(jīng)典的數(shù)值方法，結合BP神經(jīng)網(wǎng)絡得到一種新的求解雙曲守恒律方程的數(shù)值方法。數(shù)值實驗表明，該方法經(jīng)過多次訓練與網(wǎng)絡仿真，結果可自動識別斜率較大位置和間斷位置，既適用于連續(xù)問題也適用于間斷問題的求解。研究表明，數(shù)據(jù)集的選取對結果有較大影響，應選自穩(wěn)定的格式，而本文所選的熵穩(wěn)定格式在連續(xù)位置數(shù)值結果良好。該方法可根據(jù)所選多個時間層的數(shù)據(jù)預測其他待求時刻節(jié)點處的數(shù)據(jù)。