最近才入门深度学习不久,在看了arxiv.org上1902这篇文章使用神经网络在不需要已知解析解情况下就能求解常微分方程及偏微分方程的数值解,精度也很不错,自己也尝试了下,最终成功复现论文作者的结果,将代码展示一下,供需要的同学使用,才疏学浅,其中可能存在的谬误还请及时评论。
论文作者采用了一个很简单的网络结构,即只有一个单隐层的前馈神经网络(NN)且只有10个神经元,其数学表达式
m是输出的个数,这里由于是函数所以m为1,g为**函数sigmoid,w为权重,b为偏置但是这里经过测试后发现最外面的偏置项不能有,会影响最终求解的精度,所以整个网络结构就如下
是不是特别简单,这里为方便起见只画了3个神经元,但如何用这个网络结构去对我们的常微分和偏微分方程求解呢,一般来说常微分方程是如下形式
以往常规方法是先假定解的形式,如
A(x)满足边界条件,但论文作者提出了一个新观点,即不需要先假定解的形式,直接假定我们神经网络的输出就是所要求的解
而在损失函数里面加上边界条件加以限制,这样做明显简单了许多,
这里损失函数的定义就是直接将常微分方程F(x)=0左半部分直接照着写进去,对区域内每个点都计算一次F(x)并且平方之后求平均,显然如果Loss函数足够小则证明在区域上的每个点计算出来的F(x)都很趋近于0,都满足常微分方程,再加上边界条件,如一阶常微分边界条件如果是y(0)=1,则用x=0处网络的预测值y1减去1的平方加入到损失函数里面。
总结一下,损失函数里的xi就是我们在一段区域内取的点,y就是神经网络的输出,至于常微分方程中一阶以及二阶导数dy/dx如何表示,在原文中并未提到,但在其引用的[31]号文献中很清楚写明了,
k代表是几阶导数,j是输入的第j个x,i是第i个神经元,v是**函数到输出的权重,w是输入到神经元的权重,w上标k代表w的k次方,sigma代表经过**函数sigmoid的矩阵,其上方k代表sigmoid函数的k阶导数,sigmoid函数形式是,其导数很容易就能求到,令其一阶导数为y(1-y),二阶导数为y(1-y)(1-2y),所以知道这些,我们就能够将常微分方程代入求解了。我们尝试复现原论文中一个常微分方程
代码实现
import tensorflow as tf
import matplotlib.pyplot as plt
import numpy as np
x_train = np.linspace(0,2,10,endpoint=True)#生成[0,2]区间100个点
y_trail = np.exp(-0.5*x_train**2)/(1+x_train+x_train**3)+x_train**2#已知解析解用于比较
x_t = np.zeros((len(x_train),1))
for i in range(len(x_train)):
x_t[i] = x_train[i]
x1 = tf.placeholder("float", [None,1 ])#一次传入100个点[100,1]
W = tf.Variable(tf.zeros([1, 10]))
b = tf.Variable(tf.zeros([10]))
y1 = tf.nn.sigmoid(tf.matmul(x1, W)+b)#sigmoid**函数y1的形状[100,10]
W1 = tf.Variable(tf.zeros([10, 1]))
b1 = tf.Variable(tf.zeros([1]))
y = tf.matmul(y1, W1)#网络的输出[100,1]
lq = (1+3*(x1**2))/(1+x1+x1**3)
dif = tf.matmul(tf.multiply(y1*(1-y1),W),W1)#dy/dx,dif形状[100,1],即对应点的导数值
t_loss = (dif+(x1+lq)*y-x1**3-2*x1-lq*x1*x1)**2#常微分方程F的平方
loss = tf.reduce_mean(t_loss)+(y[0]-1)**2#每点F平方求和后取平均再加上边界条件
train_step = tf.train.AdamOptimizer(1e-2).minimize(loss)#Adam优化器训练网络参数
init = tf.global_variables_initializer()
sess = tf.InteractiveSession()
sess.run(init)
for i in range(50000):#训练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}))
saver = tf.train.Saver(max_to_keep=1)#保存模型,训练一次后可以将训练过程注释掉
saver.save(sess,'ckpt/nn.ckpt',global_step=50000)
saver = tf.train.Saver(max_to_keep=1)
model_file="ckpt/nn.ckpt-50000"
saver.restore(sess, model_file)
output = sess.run(y,feed_dict={x1:x_t})
output1 = sess.run(t_loss,feed_dict={x1:x_t})
y_output = x_train.copy()
y_output1 = x_train.copy()
for i in range(len(x_train)):
y_output[i] = output[i]
y_output1[i] = output1[i]
fig = plt.figure("预测曲线与实际曲线")
plt.plot(x_train,y_trail)
plt.plot(x_train,y_output)
fig2 = plt.figure("y_-y")#画实际值与预测值得偏差
plt.plot(x_train,y_trail-y_output)
fig3 = plt.figure("loss")#画出每一点对Loss的贡献
plt.plot(x_train,y_output1+(y_output[0]-1)**2)
plt.show()
结果
可以看到黄色的线是数值模拟解出的结果将蓝色的解析解的线几乎完全覆盖,说明精度很高。
这是真实值与预测值得差异,不超过0.01,所以除了一般数值模拟解微分方程的方法外,我们又多了一种利用神经网络求解的方法,增加神经元的个数和层数应该能够进一步提升求解的精度。
参考文献
https://arxiv.org/abs/1902.05563
https://arxiv.org/abs/physics/9705023