基于深度神经网络的复杂区域偏微分方程求解

郭晓斌, 袁冬芳, 曹富军*

(1. 内蒙古科技大学 信息工程学院, 内蒙古 包头 014010; 2. 内蒙古科技大学 理学院, 内蒙古 包头 014010)

偏微分方程在流体力学、气候现象、新能源、材料和航天航空领域都有着广泛的应用[1-2],研究求解偏微分方程的高精度、快速算法具有重要的意义.大多数偏微分方程都不能找到解析解,因此需要研究数值近似方法来进行求解.传统的数值方法包括有限元法、有限差分法、或有限体积法等[3-5],这些方法需要将复杂计算区域离散为计算网格,在离散网格上构造数值格式,通过求解数值格式所形成的代数方程组获得方程的近似解.尽管传统数值方法在过去的几十年中已经取得了显著的进步,并且能够处理相当复杂和高度振荡的问题,然而数值解的精度和效率严重依赖网格的质量.对于复杂计算区域问题,生成高质量的计算网格所花费的时间,在求解时间中占了相当大的比重.此外,因为高维问题中网格点的数量相对于问题维数的增加呈指数级增长,极大增加了问题求解的空间和时间复杂度,不可避免地面临“维数灾难”等问题.面对以上的挑战, 本文的目标是基于神经网络的方法求解任意形状几何域上定义的椭圆型偏微分方程. 此方法无需生成计算网格, 只需要在计算域内部和边界上进行随机采样并作为训练点集合, 因此可以处理任意形状的计算区域,避免了网格生成的困难. 此外,该方法不会受到维数增加的影响,当维数增加时,只会增加少量的参数,能够有效避免“维数灾难”问题,因此具有十分重要的研究意义和价值.

过去几十年见证了深度学习的革命,人工神经网络(ANN)的发展和进化是关键要素之一.基于人工神经网络的机器学习方法能够解决图像处理、语音识别、医学诊断等领域一些非常复杂的应用问题. 此外,人工神经网络在数学领域也已经得到了广泛的应用,因为它们能够有效地近似任意函数[6].虽然人工神经网络在几个重要的应用领域取得了令人印象深刻的成果,但关于它们为什么以及怎样有效地工作,仍然存在着许多问题.从传统意义上来讲,使用单一的隐藏层神经网络可以求解任意的偏微分方程,因为只要一个隐藏层含有足够多的神经元,就可以逼近任何函数,而且神经网络解的偏导数可以用解析闭合形式计算[7-8].研究表明,随着网络深度的增加,神经网络逼近函数的能力也越强.因此,应用深度神经网络(DNN)求解偏微分方程得到众多研究者的广泛关注.E等[9]研究了一种求解高维抛物型偏微分方程和倒向随机微分方程的新算法, 损失函数由规定的终端条件和BSDE解之间的误差给出. E等[10]提出了一种基于深度学习的深度Ritz法,通过变分的方法求解偏微分方程问题.Khoo等[11]提出了一种基于深度神经网络的求解非齐次偏微分方程的方法,将求解的物理量参数化,通过神经网络执行时间演化来寻找偏微分方程的解,进而证明使用神经网络来表征这种量的合理性.Nabian等[12]针对高维随机偏微分方程提出一种侵入性的、完全不受监督的和没有网络的基于深度学习的解决方案.Perdikaris等[13]提出了物理信息神经网络,基于训练过程实现变分问题的优化.Raissi等[14]和Owhadi[15]在机器学习法的基础上采用高斯过程拟合线性算子,并进一步将该方法推广到非线性算子的回归.Berg等[16]基于PINN的方法来求解复杂域上的偏微分方程,并取得了一些良好的结果.尽管使用PINNs在解决偏微分方程问题上有许多优点,但是PINNs仍然存在一些问题.首先,没有理论基础来了解神经网络结构的规模和所需的数据量,那么就不能保证算法不会收敛到局部最小值,并且他们的训练时间要比传统的数值方法慢.Dwivedi等[17]开发了物理信息极限学习机(PIELM),可以应用于和时间相关的线性偏微分方程问题上,同时证明了PIELM在一系列问题上超过了PINNs的准确性.Anitescu等[18]提出了一种利用人工神经网络和自适应配置策略求解偏微分方程的方法.

本文研究一种基于深度神经网络的通用偏微分方程求解方法,结合物理信息的神经网络及小样本学习求解具有复杂计算区域的问题.该方法不需要生成计算网格或设计特殊的基函数或数值格式来计算椭圆型偏微分方程的解,具有易于实现和无网格的优点,有助于求解具有复杂区域的椭圆型偏微分方程.通过对四类典型的椭圆型偏微分方程进行求解说明了该方法的效率和精度.

首先定义一个单隐层神经网络,给定d维行向量x∈Rd作为模型的输入,则神经网络的k维输出为以下形式

y=σ(xw1+b1)w2+b2

(1)

式中:w1和w2是大小为d×q和q×k的权重矩阵;b1和b2是大小为和1×k的1×q偏差向量;σ(·)为激活函数.在深度神经网络中,每一隐层通过权重矩阵和偏差对输入变量进行线性处理,然后用激活函数进行非线性变换. 上一隐层的非线性输出继续作为下一个隐藏层的输入再次重复权重矩阵的线性处理及激活函数的非线性变换.每一隐层中新的权重矩阵和偏差构成了神经网络中的一个新的隐藏层.通常,神经网络逼近复杂非线性函数的能力可以通过增加更多的隐藏层或增加隐藏层的维数来提升.具有k-1个隐藏层的神经网络可以表示为

Nθ(z)=Tk∘σ∘Tk-1∘…∘T1∘σ∘T0(z)

(2)

式中:Tk(zk-1)=wkzk-1+bk,wk表示第k-1层的与第k层之间的权重矩阵,zk-1和bk分别表示第k-1层隐藏层的输出和偏差向量,Nθ(z)为神经网络的输出,z0=z表示输入参数.σ(·)为非线性激活函数.目前,最流行的激活函数包括Sigmoid、Tanh和ReLU.ReLU激活函数是使用最广泛的函数之一,其形式为f(z)=max(0,z).但是,ReLU的高阶导数为0,这就限制了本文处理高阶导数组成的微分方程的适用性.Tanh或Sigmoid激活函数可用于二阶或高阶偏微分方程.Sigmoid激活函数是非对称的,将每个神经元的输出限制在区间[0,1]内,因此会对神经元的输出引入系统偏差.而Tanh激活函数是反对称的,通过允许每个神经元的输出在区间[-1,1]上取值,克服了Sigmoid激活引起的系统偏差问题.此外,与非对称激活的深度神经网络的训练相比,反对称激活的深度神经网络训练过程收敛的更快.

在给定多个训练数据点的回归问题中,本文可以使用欧基里德损失函数来校准权重矩阵和偏差,如下所示:

(3)

模型参数可根据式(4)进行求解

(4)

使用梯度下降法对模型进行优化,具体来说,在第i次迭代中,模型参数θ={w1,w2,…,b1,b2,…}通过式(5)更新

θ(i+1)=θ(i)-η(i)∇θJ(i)(θ(i);x,y)

(5)

式中

η(i)是第i次迭代的步长.

在反向传播过程中,利用链式法则,首先计算最后一层的梯度,最后计算第一层的梯度,当前层的部分梯度在此前的梯度计算中被重复使用.这种信息的反向流动促进了深层神经网络中每一层梯度的有效计算.算法1是梯度下降迭代算法.该算法主要用于更新神经网络中的参数.

Algorithm 1：梯度下降迭代算法

1) 给定学习率η;迭代收敛准则eps.

3) while没有达到停止准则do

采集M个样本的训练集{x1,…,xM},对应目标为yi;

计算梯度g=∇θJ(θ;x,y);

计算更新:Δθ=-η*g;

应用更新:θ←θ+Δθ;

endwhile

4) 得到神经网络训练参数θ.

根据以下微分方程求解DNN近似解u(x,θ)

(6)

式中:θ为DNN的可训练参数;N(·)是一个通用的微分算子,它是由空间导数以及线性项和非线性项组成;x是有界连续空间域D⊆RD上定义的位置向量,D中包含边界∂D.B(·)表示边界条件,由微分、线性项或非线性项组成.为了计算函数解,即计算参数θ,本文通过定义在整个离散域上的非负残差来计算,

(7)

式中:rB(θ)=rdir(θ)+rneu(θ),rdir(θ)是Dirichlet边界上的均方误差;rneu(θ)是Neumman边界上的均方误差;Nint和Nbnd分别是计算区域内部点和边界点的个数;Ndir是Dirichlet边界点的个数;Nneu是Neumman边界点的个数,且Nbnd=Ndir+Nneu.

最佳参数θ*可以根据式(8)计算,

θ*=argminr(θ)

(8)

因此, 求解微分方程(6)的解便转化为求解式(8)的优化问题,其中边界条件可以视为约束条件.根据约束条件的处理方法和严格程度,将其分为软约束和硬约束两种不同的方法.在软约束中,约束被转化为损失函数的附加惩罚项[19],这种方法易于实现,但不清楚如何调整损失函数中不同项的相对重要性,不能保证最终的解决方案会满足约束条件.在硬约束中,通过构造函数,使得任何具有该函数形式的解都保证满足边界条件[7].对于不规则的边界,构造解函数的难度很大.本文将采用以上两种约束方法计算神经网络近似解.

在软约束中,令us(x,θ)=uDNN(x,θ)为深度神经网络的近似解,根据式(9)求解DNN参数θ,

(9)

式中:λ1为惩罚项系数,用于调整损失函数中每一项的相对重要性[20].

神经网络通过式(10)更新第i次迭代中的参数

(10)

(11)

在硬约束中,近似解的函数形式可以采用以下一般形式[7]:

uh(x,θ)=C(x)+G(x,uDNN(x,θ))

(12)

式中：函数C(x)满足边界条件,函数G(x,uDNN(x,θ))不受边界条件的影响.根据式(13)求解DNN参数θ,

(13)

为了解决无约束优化问题(9)和式(13),本文采用Adam梯度下降算法.在Adam梯度下降算法的每次迭代中,损失函数的梯度使用输入空间中的一个点来近似求解,基于该点更新神经网络的参数.这种迭代更新被证明满足梯度的无偏估计[22].

同理,在硬约束中,DNN参数根据式(14)进行更新,

(14)

式中

(15)

实践中每次迭代时,损失函数的梯度是在n个不同的样本输入点上计算平均值,而不是仅在一个点上评估.

算法2总结了神经网络解法的过程.该算法可以基于预先指定的最大迭代次数停止,也可以通过损失值的动态变化作为停止标准.

Algorithm 2：复杂区域椭圆型偏微分方程的深度学习算法

输入：区域内部点:xΩ=(x1,…,xi,…,xNint)T,区域边界点:xΓ=(x1,…,xj,…,xNbnd)T.

输出：DNN中的参数值θ(i).

1) 选择求解区域,设置计算区域内部点和边界点的个数.

2) 设置DNN架构(层数,每层的维数,激活函数,以及神经网络的结构).

3) 初始化DNN参数θ(0).

4) 设置Adam优化器迭代次数Niter;设置学习率η.

5) forifrom 1 toNiterdo

(1) 根据式(2)计算出神经网络输出y;

(3) 通过式(4)和式(5)计算DNN最佳参数θ(i).

end for

本文求解复杂区域的椭圆型偏微分方程的解,大致步骤为：首先确定计算边界,并在边界上选取N个训练点.同时在计算区域内部任意选取M个训练点,将两部分训练点(M+N)作为训练集放到tensorflow框架中对神经网络模型学习,通过训练得到神经网络最优参数以及模型的近似解.然后,在计算域内选取足够的点作为测试集,并对训练后的近似解进行测试,验证该方法的准确性.

3.1 泊松方程

首先考虑具有混合边界条件的泊松方程,

其中,x,y∈[0,1].取计算区域Ω为正方形,在Ω内随机选择M=500个训练点,在∂Ω上均匀选择N=200个训练点,将M和N共同作为训练集,如图1所示.且方程的解析解为

图1 训练集Fig.1 The traning set

u(x,y)=e-x(x+y3)

损失函数构造如下:

其中,λ1为1 000,λ2为30,选择学习率η为0.001,经过30 000次的迭代训练,得到方程的训练解.取M=200 000,N=2 000作为测试集,测试集的解析解和DNN近似解如图2、图3所示.

图2 解析解Fig.2 The analysis solution

图3 DNN近似解Fig.3 The DNN solution

图4给出了测试集上的误差,可以看出神经网络求解的结果与解析解的误差在-1%～5‰之间,效果比较理想.

图4 测试集误差Fig.4 The error of the testing sets

针对泊松方程,分别采用Adam,SGD,Adadelta,RMSProp,Adagrad等5种梯度优化算法处理损失函数,损失函数收敛曲线如图5所示.横轴代表的是梯度下降算法迭代次数,纵轴表示损失函数的损失值.随着迭代次数的增加,损失值逐渐减小.其中,Adam优化算法对模型的优化效果要比其余方法要好,主要是因为该算法结合了Adagrad和RMSProp算法的优点,相对于传统的梯度下降算法来说,Adam优化算法具有自适应能力,计算收敛更快.SGD优化算法结果与Adam算法的梯度下降曲线效果相近,但是,当数据集比较大的时候,SGD算法的收敛速度和效率优于Adam,所以在不同的情况下应该选用合适的优化算法.本文计算表明,采用Adam优化算法可以达到理想的效果.

图5 不同优化方法的收敛性Fig.5 The convergence of different optimization methods

3.2 二维线形平流方程

平流方程在大气运动的研究中具有重要意义,是大气运动中重要的过程之一,运动方程、热通量方程和水汽方程中都含有平流项.然而,数值求解方法很难求出此类方程的解析解,因此研究神经网络解平流方程具有非常重要的现实意义.

其中:x,y∈[0,1],且满足混合边界条件,方程的解析解为

u(x,y)=sin(3πx)cos(3πy)

其中,a,b是平流系数,在此问题中,取a=2,b=3,取计算区域Ω为正方形,在Ω内生成M=500个任意的训练点,在∂Ω上生成N=200个训练点,如图6所示.

图6 训练集Fig.6 The traning set

图7、图8分别为测试集在解析解和DNN解上的函数结果,图9给出了二维线性平流方程在测试集上的误差,误差结果在-4‰～8‰之间.

图7 解析解Fig.7 The analysis solution

图8 DNN近似解Fig.8 The DNN solution

图9 测试集误差Fig.9 The error of the testing set

损失函数构造如下:

其中:λ1=200,λ2=50,λ3=50,Nneu1是y=1时的Neumman边界点的个数,Nneu2是y=0时的Neumman边界点的个数.取M=200 000,N=2 000作为测试集.

图10比较了随着迭代次数的增加,不同优化算法的收敛曲线.

图10 不同优化方法的收敛性Fig.10 The convergence of different optimization methods

3.3 二维对流扩散方程

本算例考虑二维对流扩散方程

其中:计算区域为花瓣形,花瓣的半径[23]可以表示为Φ(x,y)=r-(0.4+0.3sin(5θ)),其中

将Ω区间划分为80×80个网格点,判断网格点是否落在花瓣半径以内,当-0.1<Φ(xi,yi)<0时,此时的(xi,yi)i=1,2,…,M是花瓣的边界,当Φ(xj,yj)<-0.1,(xj,yj)j=1,2,…,N为花瓣的内部点,如图11所示.

图11 训练集Fig.11 The traning set

方程的解析解为

损失函数构造如下:

其中:λ1为1 000选取花瓣内的网格点M=14 400,N=6 300作为测试集.图12和图13给出了测试集在解析解和DNN近似解函数图像.图14给出了测试集的误差,误差范围在-6‰～4‰之间,表现出较高的精度.图15比较了随着迭代次数的增加,不同优化方法的损失函数收敛曲线,可以看出Adam方法的收敛效率优于其它四种方法.

图12 解析解Fig.12 The analysis solution

图13 DNN近似解Fig.13 The DNN solution

图14 测试集误差Fig.14 The error of the testing set

图15 不同优化方法的收敛性Fig.15 The convergence of different optimization methods

3.4 二维势能方程——由梯度势能计算力场

在分子动力学中,力场可以由势能的负梯度求得

方程的解析解为

u(x,y)=1+y2-2x2+x6

损失函数构造如下：

其中：λ1为100,F(x,y)是力场,u(x,y)是势能.计算区域Ω为中国大陆,在Ω内任意选择M=500个训练点,在∂Ω上选择N=455个训练点,如图16所示.采用具有2个隐藏层、30个神经元的神经网络计算训练解.取M=100 000,N=2 276作为测试集代入ANN训练解中.图17为测试集上的精确解.图18为测试集上的近似解.图19给出了测试集上的误差,误差范围在-1%～7‰之间.图20比较了随着迭代次数的增加,不同优化方法的损失收敛曲线.由于不同偏微分方程的复杂度是不同的,复杂度高的偏微分方程在训练过程中不容易训练,损失函数收敛过程比较慢,且不易收敛,在训练过程中需要增加训练次数,并采用适当的优化方法来处理.

图16 训练集Fig.16 The traning set

图17 解析解Fig.17 The analysis solution

图18 DNN近似解Fig.18 The DNN solution

图19 测试集误差Fig.19 The error of the testing set

图20 不同优化方法的收敛性Fig.20 The convergence of different optimization methods

本文研究一种利用神经网络求解复杂区域椭圆型偏微分方程的方法.该方法使用完全无监督的类梯度下降算法进行迭代训练,在求解过程中无需调用其他的偏微分方程数值解算器.此外,它在整个空间域上是无网格的,因此,有利于处理不规则的计算区域.通过反向传播算法的推导,将二阶偏微分方程扩展到任意阶偏微分方程组.数值算例结果表明,在本文提供的神经网络计算框架下,任意椭圆型偏微分方程在任何区域内都可以通过训练得到它们的近似解,且具有较高的精度;同时,神经网络方法求解复杂区域椭圆型偏微分方程的收敛性和计算效率还需要进一步提升和改进.因此,还需要考察不同的优化方法,进一步提高神经网络方法的性能.

致谢：本文得到内蒙古科技大学创新基金(2019YQL02)项目的资助,在此表示感谢.

猜你喜欢梯度方程网格方程的再认识中学生数理化·七年级数学人教版(2022年5期)2022-06-05一个带重启步的改进PRP型谱共轭梯度法数学物理学报(2022年1期)2022-03-16一个改进的WYL型三项共轭梯度法数学物理学报(2021年6期)2021-12-21方程(组)的由来中学生数理化·七年级数学人教版(2021年5期)2021-11-22随机加速梯度算法的回归学习收敛速度数学物理学报(2021年5期)2021-11-19圆的方程新世纪智能(数学备考)(2020年12期)2020-03-29追逐作文新天地(初中版)(2019年6期)2019-08-15一个具梯度项的p-Laplace 方程弱解的存在性华东师范大学学报（自然科学版）(2019年3期)2019-06-24重叠网格装配中的一种改进ADT搜索方法北京航空航天大学学报(2017年6期)2017-11-23基于曲面展开的自由曲面网格划分浙江大学学报（工学版）(2016年10期)2016-06-05