楊 震,陳豫眉,李 霜

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

0 引言

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

1 預(yù)備知識(shí)

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

圖1 人工神經(jīng)元的結(jié)構(gòu)模型Fig.1 Structural Model of Artificial Neuron

上述神經(jīng)元工作過程的數(shù)學(xué)模型如下：

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

階躍函數(shù)將神經(jīng)元的輸入映射為0或1,分別代表神經(jīng)元的“抑制”與“興奮”兩種狀態(tài).但是,由于階躍函數(shù)在階躍點(diǎn)處的不連續(xù)和不可導(dǎo),實(shí)際應(yīng)用中常常采用圖2(b)所示的Sigmoid函數(shù)：

Sigmoid作為激活函數(shù),它能將神經(jīng)元的輸入值連續(xù)的映射到(0,1)范圍內(nèi),并且能夠進(jìn)行多次求導(dǎo),符合概率分布的特點(diǎn).

2 單隱層神經(jīng)網(wǎng)絡(luò)的網(wǎng)絡(luò)結(jié)構(gòu)與BP算法

圖3 單隱層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.3 Network Structure of Single Hidden Layer Neuron

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

圖4 節(jié)點(diǎn)j到輸出層的信號(hào)流圖Fig.4 Signal Flow Diagram of Node j to Output Layer

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

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

其中,η為學(xué)習(xí)率. 第n次迭代結(jié)束后,輸出層權(quán)值Wm+j(n)調(diào)整為：

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

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

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

根據(jù)梯度下降法可知,權(quán)值Wj和閾值bj的改變量分別為：

其中,η為學(xué)習(xí)率. 第n次迭代結(jié)束后,輸出層權(quán)值Wj和閾值bj分別調(diào)整為：

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

圖5 輸出層經(jīng)過節(jié)點(diǎn)j到輸出層的信號(hào)流圖Fig.5 Signal Flow Diagram of the Output Layer Passing Through the Node j to the Output

3 神經(jīng)網(wǎng)絡(luò)求解常微分方程的算法思路

在實(shí)際問題中,期望值或真實(shí)值往往未告知,對(duì)于人工神經(jīng)網(wǎng)絡(luò)誤差e的定義需要作出相應(yīng)的調(diào)整.一般的n階常微分方程具有形式

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

定義損失函數(shù)Loss為

4 網(wǎng)絡(luò)推演過程及python程序

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

采用python語言求解常微分方程的近似解,其優(yōu)點(diǎn)是處理大數(shù)據(jù)更加方便快捷,運(yùn)行速度較快.程序編寫可采用以下五個(gè)步驟實(shí)現(xiàn)：

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

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

(3)以損失函數(shù)是否滿足停機(jī)條件來判斷神經(jīng)網(wǎng)絡(luò)權(quán)值和閾值的更新；若不滿足,則通過反向傳播依次更新權(quán)值和閾值；

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)經(jīng)過一次反向傳播,更新權(quán)值和閾值,然后繼續(xù)重復(fù)第(2)～(3)步；

(5)最終得出滿足停機(jī)條件的神經(jīng)網(wǎng)絡(luò).

4 數(shù)值算例

表1 幾種方法對(duì)應(yīng)的計(jì)算結(jié)果Tab.1 The Corresponding Results of Several Methods

表2 幾種方法的計(jì)算結(jié)果Tab.2 The Results of Several Methods

通過數(shù)值算例,本文介紹的單隱層神經(jīng)網(wǎng)絡(luò)解法對(duì)于常微分方程問題的計(jì)算結(jié)果與改進(jìn)歐拉方法,Heun's方法和文獻(xiàn)[3]中方法相比,計(jì)算精度更高,誤差更小.并且誤差可以通過增加隱層神經(jīng)元的個(gè)數(shù)或?qū)訑?shù),或者使用遺傳算法優(yōu)化權(quán)重等方法去縮小,因此常微分方程可用神經(jīng)網(wǎng)絡(luò)求解.

5 結(jié)語

本文討論了一種單隱層神經(jīng)網(wǎng)絡(luò)算法在數(shù)值常微分方程求解中的應(yīng)用.文中將微分方程的求解問題與損失函數(shù)的優(yōu)化問題聯(lián)系起來,以損失函數(shù)的優(yōu)化問題去訓(xùn)練神經(jīng)網(wǎng)絡(luò)進(jìn)而得出微分方程的近似解,并通過網(wǎng)絡(luò)推演觀察到神經(jīng)網(wǎng)絡(luò)權(quán)值跟閾值的變化情況,最后使用python語言求解的數(shù)值算例證明了神經(jīng)網(wǎng)絡(luò)算法的精確性和實(shí)用性.