基于python語言的一種常微分方程神經網絡解法

楊 震,陳豫眉,李 霜

(1.西華師范大學數學與信息學院，四川南充 637009；2.西華師范大學公共數學學院，四川南充 637009；3.西華師范大學計算方法及應用軟件研究所,四川南充 637009)

0 引言

在科學計算及其工程應用中,常常需要求解微分方程.但解析方法只能用來求解某些特殊類型的方程,其他情況主要依靠數值解法.對于常微分方程問題,常用的差分方法有Euler方法,改進Euler方法和霍因(Heun)方法等[1，2].其中,Euler方法是一種最簡單的單步法,計算量小,但精度低.改進的Euler方法由于增加了校正過程,在解的精度上有了較大的改進. Heun's方法是在Euler方法的基礎上進一步改進,精度更高.近年來,神經網絡在機器學習,深度學習,人工智能等各個領域取得了顯著的發展,應用神經網絡解決微分方程的問題已經成為熱門.文獻[3]中利用余弦基神經網絡求解常微分方程.文獻[4]中利用神經網絡與小波分析求解微分方程.1988年,Lagaris和Likas提出了一種求解偏微分方程初邊值問題的神經網絡算法[5].由于單隱層神經網絡可按任意精度去逼近任何的連續函數,所以該算法可以用于解決微分方程問題.2010年,Baymani,Kerayechia,EffatiS使用神經網絡算法去求解stokes方程[6].2017年,Raissi,Perdikaris,Karniadakis提出基于高斯過程的神經網絡算法[7].Long,Lu和Ma提出基于卷積微分算子求解偏微分方程的神經網絡算法[8].2018年,Han,Jentzen等人提出了一種使用多層神經網絡求解高維向后微分方程的算法[9].2019年,He和Xu討論了卷積神經網絡與多重網格之間的關系[10].本文所采用的神經網絡算法將常微分方程的殘差及初邊值條件殘差之和作為損失函數,通過最小化損失函數求解和評估模型來訓練神經網絡中的權重,從而得出需要點的近似解.基于python語言,Adam優化器來解決離散優化問題.

1 預備知識

1943年,Mcculloch,Pitts將生物神經元模型抽象為圖1所示的人工神經元模型[11].設神經元的輸入為x1,x2,…,xn,這些輸入通過帶有權值w1,w2,…,wn的連接來進行傳遞, 然后通過加權求和得到總輸入,再與神經元的閾值b作出比較,最后利用激活函數σ進行非線性變化得到神經元的輸出y.人工神經元結構模型如圖1所示：

圖1 人工神經元的結構模型Fig.1 Structural Model of Artificial Neuron

上述神經元工作過程的數學模型如下：

圖2(a) 階躍函數圖像Fig.2(a) Graph of Step Function圖2(b) Sigmoid函數圖像Fig.2(b) Graph of Sigmoid Function

階躍函數將神經元的輸入映射為0或1,分別代表神經元的“抑制”與“興奮”兩種狀態.但是,由于階躍函數在階躍點處的不連續和不可導,實際應用中常常采用圖2(b)所示的Sigmoid函數：

Sigmoid作為激活函數,它能將神經元的輸入值連續的映射到(0,1)范圍內,并且能夠進行多次求導,符合概率分布的特點.

2 單隱層神經網絡的網絡結構與BP算法

圖3 單隱層神經網絡結構Fig.3 Network Structure of Single Hidden Layer Neuron

神經網絡的學習算法有很多,針對本文的單隱層神經網絡問題,最常用的就是BP算法,(即誤差反向傳播算法)[12].BP算法由信號的正向傳播與誤差的反向傳播兩個部分組成.信號的正向傳播即輸入值通過輸入層進入神經網絡,經過隱藏層逐步傳遞至輸出層,如果輸出值跟期望值不符,則進入誤差的反向傳播過程；如果輸出值跟期望值相符,則結束學習算法.誤差的反向傳播即將誤差信號(輸出值與期望值之差)按原來的連接反傳計算,經過隱藏層最終達到輸入層,在反向傳播的過程中將誤差分攤給每層的每個單元,得到每層每個單元的誤差信號,并且將其作為修正每個單元權重的依據.這個過程一般采用梯度下降法[13]來完成,在不停修正各神經元的權值跟閾值的過程中,誤差信號逐漸減小.

圖4 節點j到輸出層的信號流圖Fig.4 Signal Flow Diagram of Node j to Output Layer

則e(n)對Wm+j(n)的梯度為：

根據梯度下降法可知,權值Wm+j(n)的改變量為：

其中,η為學習率. 第n次迭代結束后,輸出層權值Wm+j(n)調整為：

Wm+j(n+1)=Wm+j(n)+ΔWm+j(n).

從而可得e(n)對Wj(n)的梯度為：

則e(n)對bj(n)的梯度為：

根據梯度下降法可知,權值Wj和閾值bj的改變量分別為：

其中,η為學習率. 第n次迭代結束后,輸出層權值Wj和閾值bj分別調整為：

Wj(n+1)=Wj(n)+ΔWj(n),bj(n+1)=bj(n)+Δbj(n).

圖5 輸出層經過節點j到輸出層的信號流圖Fig.5 Signal Flow Diagram of the Output Layer Passing Through the Node j to the Output

3 神經網絡求解常微分方程的算法思路

在實際問題中,期望值或真實值往往未告知,對于人工神經網絡誤差e的定義需要作出相應的調整.一般的n階常微分方程具有形式

考慮如下的常微分方程初邊值問題：求解y,滿足

定義損失函數Loss為

4 網絡推演過程及python程序

圖6 第n次推演過程圖Fig.6 Diagram of the nth Derivation

采用python語言求解常微分方程的近似解,其優點是處理大數據更加方便快捷,運行速度較快.程序編寫可采用以下五個步驟實現：

(1)以包含邊值條件的損失函數Loss滿足收斂精度要求或給定迭代次數為停機條件,隨機選取初始權值和閾值；

y1 = tf.nn.sigmoid(tf.matmul(x1, W)+b) W1 = tf.Variable(tf.zeros([g, 1]))y = tf.matmul(y1, W1)

(3)以損失函數是否滿足停機條件來判斷神經網絡權值和閾值的更新；若不滿足,則通過反向傳播依次更新權值和閾值；

train_step = tf.train.AdamOptimizer(1e-2).minimize(Loss) init = tf.global_variables_initializer()sess = tf.InteractiveSession()sess.run(init)for i in range(50000): sess.run(train_step, feed_dict={x1: x_t}) if i % 50 == 0: total_Loss = sess.run(Loss, feed_dict={x1: x_t}) print("Loss={}".format(total_Loss)) print(sess.run(y[0], feed_dict={x1: x_t}))

(4)經過一次反向傳播,更新權值和閾值,然后繼續重復第(2)～(3)步；

(5)最終得出滿足停機條件的神經網絡.

4 數值算例

表1 幾種方法對應的計算結果Tab.1 The Corresponding Results of Several Methods

表2 幾種方法的計算結果Tab.2 The Results of Several Methods

通過數值算例,本文介紹的單隱層神經網絡解法對于常微分方程問題的計算結果與改進歐拉方法,Heun's方法和文獻[3]中方法相比,計算精度更高,誤差更小.并且誤差可以通過增加隱層神經元的個數或層數,或者使用遺傳算法優化權重等方法去縮小,因此常微分方程可用神經網絡求解.

5 結語

本文討論了一種單隱層神經網絡算法在數值常微分方程求解中的應用.文中將微分方程的求解問題與損失函數的優化問題聯系起來,以損失函數的優化問題去訓練神經網絡進而得出微分方程的近似解,并通過網絡推演觀察到神經網絡權值跟閾值的變化情況,最后使用python語言求解的數值算例證明了神經網絡算法的精確性和實用性.