張 慊 喬 丹 唐少強,2)

*(北京大學工學院力學與工程科學系，北京100871)

?(北京大學數(shù)學科學學院概率統(tǒng)計系，北京100871)

分子動力學計算常用于精細分析材料的性質，往往采用周期邊界條件。對于多尺度計算和其他一些力學研究，周期邊界條件不適用，需要設計人工邊界條件[1]。本文將針對一維線性原子鏈，運用機器學習設計人工邊界條件。該方法需要的先驗知識很少，且可以快速產生一系列邊界條件。

傳統(tǒng)的獲得人工邊界條件的方法有時間歷史積分[2]與在給定波數(shù)附近進行泰勒展開的方法[3]等。本文提出的方法主要是受后者的啟發(fā)。

1 原子鏈的人工邊界條件

考慮一維無窮長原子鏈，每個原子與其鄰近原子線性相互作用。歸一化控制方程為

其中uj表示第j個原子的位移。

方程(1)的解可表示為模態(tài)

的線性疊加，其中

為色散關系。ξ∈[-π,0]為右行波，ξ∈[0,π]為左行波。

如果數(shù)值截斷模擬有限長度內的原子鏈，需要在數(shù)值邊界上設計人工邊界條件。以左邊界u0為例，我們設計以下形式的線性邊界條件

N表示邊界條件用到的原子數(shù)，我們取N= 6,c0= 1。其中cj及bj(j= 0,1,···,6) 即為待定的邊界條件系數(shù)。

定義殘差函數(shù)

可以發(fā)現(xiàn)如果波數(shù)為ξ的左行波能夠完全滿足邊界條件，則Δ(ξ) = 0。而如果不能完全滿足邊界條件，Δ(ξ) /= 0，就會產生反射?？梢灾婪瓷湎禂?shù)

2 前饋神經網絡

前饋神經網絡是神經網絡的一種，由一個輸入層、若干個隱藏層和一個輸出層組成。每一層由一組權重和一個激活函數(shù)組成，上一層的輸出經過權重的線性組合和激活函數(shù)的映射之后輸出給下一層隱藏層，直到最終輸出。

最簡單的前饋神經網絡結構如圖1 所示(這是沒有隱藏層的情況)：

圖1 單層神經網絡

神經網絡學習從輸入

到輸出的映射，即學習函數(shù)hW,b(x)滿足

其中θ為參數(shù)，代表W,b。

為對應每個xi的權重，b為常數(shù)項。f稱為激活函數(shù)，在得到輸入加權的結果后，用激活函數(shù)將該值映射到最終的輸出。常見的激活函數(shù)有：線性函數(shù)(相當于不使用激活函數(shù))，ReLU，sigmoid(用于二分類問題)，softmax(用于多分類問題)，Leaky ReLU，tanh等。

兩層的前饋神經網絡如圖2所示。

圖2 多層神經網絡

更一般地，輸出也可以是一個向量。

訓練神經網絡一般使用反向傳播算法，本質是一種梯度下降算法，即通過求損失函數(shù)J關于每個參數(shù)的偏導數(shù)，然后用梯度下降更新參數(shù)。損失函數(shù)J可以定義成均方誤差

梯度下降算法即α是學習率，一般設置為一個較小的正數(shù)。其中h(X;θ)表示當前參數(shù)下神經網絡的輸出，?h(X;θ)表示訓練的目標值。

對多層的前饋神經網絡，使用鏈式法則求損失函數(shù)J關于各個參數(shù)的梯度。記是第k-1 層第j個神經元連接到第k層的第i個神經元的權值表示第k層第j個神經元的偏置表示第k層第j個神經元的輸入表示第k層第j個神經元的輸出，即為激活函數(shù)。

第k層第j個神經元的誤差(實際值與預測值的誤差)定義為應用鏈式法則，有

其中

對上式兩邊微分

整體梯度下降以及隨機梯度下降的收斂速度通常較慢，可使用Adam算法[4]進行優(yōu)化。

3 數(shù)值算例

對于給定初始條件u=fi(x)，設fi(x)支集為閉區(qū)間E= [-L/2,L/2]。在區(qū)間E0= [-Lc/2,Lc/2]上用二階Verlet 算法對時間離散進行計算，將總長度取得充分長(Lc＞5L)并計算足夠步數(shù)(使波前不超出邊界)，選取計算區(qū)間中一個點作為左邊界點，記錄對應的u0,···,uN-1˙u0,···,˙uN-1，即得參考數(shù)值解，用于訓練一個單層的前饋神經網絡。該神經網絡以˙u0為輸出，以u0,···,uN-1˙u1,···,˙uN-1為輸入。由線性約束假設(4)，函數(shù)應該是線性的，因此選取激活函數(shù)為線性激活函數(shù)，即f(x)=x。因為區(qū)間上原子位移全部為零對應于原子鏈靜止的情形，邊界原子的速度也一定是零，因此零輸入一定對應零輸出，所以偏置b取為零。于是訓練得到的神經元的權值就是邊界條件(4)中的系數(shù)。

以N= 6 為例，取L= 120,Lc= 1200，令jlb=-90處作為左邊界點。取初值為

A為任意正常數(shù)，M為正整數(shù)，這里令A= 50，M= 17。離散時間步長Δt= 0.005，則單位時間有200步，計算400個單位時間。將作為訓練集的輸入作為訓練集的輸出。訓練200 輪之后得到系數(shù)cjbj,j= 0,1,2,···,6(如表1所示)。與文獻[3]提出的邊界條件(表2 所示)對比，兩者差別很大。為了驗證這組系數(shù)的正確性，在區(qū)間[-50,50]上取初值

然后分析機器學習方法得到的邊界條件的殘差。計算Δ關于ξ的關系，結果如圖4所示。

從圖4 中可見，文獻[3]提出的人工邊界條件的殘差在ξ→π，即短波處發(fā)散到無窮。而機器學習得到的人工邊界條件雖然在ξ→0，即長波處殘差較大，但是短波處殘差是有界的。

為進一步考察機器學習得到的人工邊界條件在能量意義下的精確度，計算邊界條件對應的反射系數(shù)R關于ξ的關系，結果如圖5 所示。這里有反射系數(shù)大于1 的波段，但前述算例仍是穩(wěn)定的；此外，還可以通過在訓練數(shù)據(jù)中增加更多這個波段的模態(tài)來改善。

表1 N =6 時機器學習得到的人工邊界條件的系數(shù)

表2 N =6 時文獻[3]提出的人工邊界條件的系數(shù)

圖3 N =6 時機器學習方法與參考解的數(shù)值結果比較，實線代表參考解，圓圈代表采用人工邊界條件的數(shù)值解。由于解的對稱性，只畫出了通過右側邊界的波形

圖4 N =6 時Δ 關于ξ 的關系，其中點線與ξ 軸重合，表示精確的邊界條件，點劃線代表文獻[3]提出的人工邊界條件的殘差，實線代表機器學習得到的人工邊界條件的殘差

圖5 N =6 時R 關于ξ 的關系，其中點線與ξ 軸重合，表示精確的邊界條件，點劃線代表文獻[3]提出的人工邊界條件的反射系數(shù)，實線代表機器學習得到的人工邊界條件的反射系數(shù)

4 結論

對于線性原子鏈的人工邊界條件這一問題，本文在基于泰勒展開的匹配邊界條件方法基礎上，提出了基于前饋神經網絡的機器學習算法。該方法不需要任何對方程(1)的細致分析就能得到一組較好的邊界條件。雖然該方法在一維線性原子鏈模型上的應用仍然存在殘差和反射系數(shù)較大的問題，但是應用該方法需要的先驗知識很少，編程方便且容易實現(xiàn)，因此具有較高的實用價值，可以進一步推廣到其他人工邊界條件的設計。