热传导界面问题的一种神经网络算法

王铸础,谢小平

(四川大学数学学院,成都 610064)

1 引 言

近年来，神经网络在机器学习、人工智能领域取得了极大进展，被广泛应用于计算机视觉，自然语言处理，强化学习等诸多方面. 其中应用神经网络算法求解偏微分方程已成为一个热门的研究领域[1-17].

1998年，Lagaris和Likas[9]提出了一种用于求解偏微分方程初边值问题的人工神经网络(ANN)算法.该算法将近似解表示为两项之和：第一项的函数满足初边值条件；第二项为一个在满足在初始及边界上恒等于0的函数与一个前馈神经网络函数的乘积，并通过使用前馈神经网络对连续函数的逼近性能，以使近似解逼近方程真解.为训练第二项中的神经网络， 该算法将其权重作为自变量，以极小化一个适当的损失函数.由于单隐层神经网络可以以任意精度逼近连续函数[2]，所以该类神经网络可被用于求解微分方程.文献[10]讨论了边界不规则时的推广算法.2010年，Baymani,Kerayechian及Effati使用文献[9]中的神经网络算法求解了Stokes方程[1].2017年，Raissi,Perdikaris,Karniadakis 提出了一种基于高斯过程的神经网络算法[15]，并用其求解线性微分方程.2017年，Long，Lu和 Ma提出了求解偏微分方程的一种所谓“Pde-net”的神经网络[12]，其特点是使用卷积核学习微分算子.2018年，Han，Jentzen和E[4]提出了一种通过使用多层神经网络计算随机游走项相关梯度，以求解高维向后随机微分方程的算法.基于Ritz 法与神经网络，E和Yu 提出了Deep Ritz 方法[3]；He，Li，Xu和Zheng 从使用RELU 激活函数的神经网络函数是分片线性多项式这一角度出发，证明了自由度为N的d维线形有限元函数可由至多[lognd]+1层ReLu神经网络表出[5].Sirignano和 Spiliopoulos提出了一种求解偏微分方程的深度学习算法，其损失函数为方程残差与初边值条件残差的L2范数的平方和[16].2019 年，He和Xu在文献[6]中讨论了卷积神经网络与多重网格之间的关系.

本文将研究Lagaris和Likas所提出的神经网络算法[9]在热传导界面问题中的推广应用. 我们在每个子区域构造形如A(x,t)+B(x,t)*N(x,t)的近似解，其中A(x,t) 满足原方程的初边值条件，B(x,t) 为一在初始时刻和边界上恒等于零的函数，N(x,t) 为神经网络函数.对于界面条件的处理，我们将Dirichlet 界面条件视为边界条件强制满足，将Neumann界面条件的残差作为损失函数的一部分.最后，通过极小化损失函数，我们得到优化的神经网络函数权重，从而最终获得问题的近似解.针对一维问题，我们使用Adam 算法求解离散优化问题.针对高维情形，我们给出了一种基于随机梯度法思想的类随机梯度法来求解相应的离散优化问题.本文的神经网络算法既保留了原算法近似解表达形式简洁、适用范围广等优点，又易于处理界面条件，便于推广应用于其它类型的界面问题.

2 预备知识

令N+为正整数集.给定n,m∈N+，定义线性函数

(1)

(2)

给定整数s,r≥1和矩阵[Qij]∈Rs×r,定义

(3)

(4)

αi,bi,β∈R,Wi∈R1×d,i=1,2,…,M}

(5)

注1为方便起见，后文用p表示式(5)中的参数αi,bi,Wi,β，并记

p=(α1,…,αM,b1,…,bM,W1,…,WM,β).

我们称其为神经网络函数的权重.

3 算 法

设Ω为有界区域，Ω⊂Rd,d=1,2,3. 一个光滑的界面将区域Γ分成Ω1,Ω2两部分.设n1,n2分别为Ω1,Ω2在界面Γ上的单位外法向量.设时间长度T>0.考虑如下的热传导初边值界面问题：求解u(x,t),v(x,t)，满足

ut-▽·(a1▽u)=f1,(x,t)∈Ω1×(0,T]

(6)

ut-▽·(a2▽u)=f2,(x,t)∈Ω2×(0,T]

(7)

u=g1,(x,t)∈∂Ω1Γ×(0,T]

(8)

u=g2,(x,t)∈∂Ω2Γ×(0,T]

(9)

u(x,0)=u0,x∈Ω1

(10)

v(x,0)=v0,x∈Ω2

(11)

u-v=φ,(x,t)∈Γ×(0,T]

(12)

a1▽u·n1+a2▽v·n2=ψ,(x,t)∈Γ×(0,T]

(13)

其中a1,a2>0为热传导系数，f1,f2为源项，f1,f2,g1.g2.u0.v0,φ,ψ均为各自定义域上的连续函数.

引入损失函数

Lossu(u(x,t))=(ut(x,t)-

▽·(a1▽u(x,t))-f1(x,t))2,

Lossv(v(x,t))=(vt(x,t)-

▽·(a2▽v(x,t))-f2(x,t))2,

LossJ(u(x,t),v(x,t))=(a1▽u(x,t)·n1+

a2▽v(x,t)·n2-ψ(x,t))2.

给定参数α>0,定义如下优化问题：

(14)

其中

D1={(x,t)|x∈Ω1,t∈[0,T]},

D2={(x,t)|x∈Ω2,t∈[0,T]},

D3={(x,t)|x∈Γ,t∈[0,T]}.

易知,原问题(式6～13)的真解为此优化问题的最优解.本文的神经网络算法的基本思路是：首先利用神经网络函数构造满足初边值条件和Dirichlet界面条件的近似解U,V,然后求解如下的关于神经网络权重的离散优化问题：

(15)

其中Xi⊂Di(i=1,2,3)为给定的离散点集(训练集).

下面给出近似解U,V的构造步骤.

第一步，设

U(x,t,p1)=A1(x,t)+

B1(x,t)*N1(x,t,p1),(x,t)∈Ω1×[0,T]

(16)

其中A1和B1满足

A1(x,t)=g1,(x,t)∈∂Ω1Γ×(0,T],

A1(x,0)=u0,x∈Ω1,

B1(x,t)=0,(x,t)∈∂Ω1Γ×(0,T],

B1(x,0)=0,x∈Ω1,

β(1),x∈Ω1,t∈[0,T]

(17)

第二步，设

V(x,t,p2)=A2(x,t)+B2(x,t)*

N2(x,t,p2),(x,t)∈Ω2×[0,T]

(18)

其中A2和B2满足

A2(x,0)=v0,x∈Ω2,

B2(x,t)=0,(x,t)∈∂Ω2×(0,T],

B2(x,0)=0,x∈Ω2,

将上述U,V代入离散优化问题(15)，则离散优化问题转化为关于权重p1,p2的优化问题，并通过优化算法进行求解.

注2构造近似解U的关键在于构造满足关于u的初边值条件的函数A1，构造近似解V的关键在于构造满足关于v的初边值条件与Dirichlet 界面条件的函数A2.

为验证算法的有效性，下面定义三种形式的误差:

其中|X1|,|X2|,|X3|分别表示优化问题(14)中离散点集X1,X2,X3包含的点的个数.

源项

f1(x,t)=0.09πe(-0.1π2t)(cos(πx)-

f2(x,t)=0.06π2e(-0.1π2t)sin(4πx),

初值条件

边值条件

g1(0,t)=g2(1,t)=0,t∈[0,t]，

界面条件

方程的真解为

u(x,t)=e(-0.1π2t)sin(πx),

v(x,t)=e(-0.1π2t)sin(4πx),

按照上一节给出的神经网络算法步骤，我们首先构造近似解

U(x,t) = (1-t)*sin(πx) +

x*t*N1(x,t,p1),

V(x,t)=(1-t)*sin(4πx)+

(1-x)*(3x-2)*t*N2(x,t,p2),

其中神经网络N1,N2的隐层神经元数M1=M2=10. 离散优化问题(14)中的离散训练点集取为：

0

0

将U,V代入式(15)，取参数α=1，然后使用Adam算法[7]求解，步长取0.01，其余可调参数均取默认值；神经网络权重初值随机生成.

表1给出了误差ER,E2和EMS随Adam 算法迭代次数变化情况，图1给出了相应的变化曲线.特别地，右图实线、虚线分别为测试集(见注4)与训练集 上的EMS误差随迭代次数变化的曲线.可以看出：随着迭代次数增加，近似解精度逐步提高；在迭代近100次后，误差趋于稳定；右图的两条误差变化曲线基本一致，表明我们的算法具有较好的泛化能力.

表1 不同误差随迭代变化情况：算例 4.1

图2至图9给出了近似解图像随迭代次数变化的对比，迭代次数的选取与表1 对应，真解图像见图11. 虽然初始误差较大，但随着迭代次数增加，近似解越来越接近真解，这表明本文的神经网络算法具有学习界面问题真解的能力.

图2 迭代初值Fig.2 The initial value

图3 迭代16次Fig.3 16 times iteration

图4 迭代 24 次Fig.4 24 times iteration

图5 迭代32次Fig.5 32times iteration

图6 迭代48次Fig.6 48 times iteration

图7 迭代 72 次Fig.7 48 timesiteration

图8 迭代108次Fig.8 108 times iteration

图9 迭代148次Fig.9 148 times iteration

图10 真解Fig.10 The exact solution

注4本算例中测试集选取为：

0

0

记测试集上的EMS误差为

为考察了离散优化问题(15)中参数α的选择对近似解精度的影响，我们分别取α=0.1,1,10,100,1 000,10 000.图11给出了相应的E2误差随迭代变化的曲线，其中因为α=0.1,1,10,对应的三条曲线基本重合，所以在图中只给出α=0.1的曲线作为这三种情形的代表.

图11 不同参数α 对应的E2 误差随迭代次数的变化曲线: 算例4.1，α=0.1,100,1 000,10 000

从图11可以看出，对所有情况，随着迭代次数增加，近似解的E2误差均逐步减小； 到达一定迭代次数后，误差趋于稳定.特别地，当α=100时，其对应的误差在迭代100次后小于其它参数选择对应的误差.综合考虑精度与速度，α=100 是较优的选择.

算例4.2考虑二维常系数问题(6)～(13)，其中Ω1=[0,0.5)×[0,1],Ω2=(0.5,1]×[0,1],T=1，热传导系数

a1(x,y)=1,(x,y)∈[0,0.5]×[0,1],

a2(x,y)=0.5,(x,y)∈(0.5,1]×[0,1],

源项

f1(x,y,t)=e-t((x-0.5)sin(x+y)-2cos(x

+y)),(x,y,t)∈[0,0.5]×[0,1]×[0,1],

f2(x,y,t)=-2e-tcos(x+y),(x,y,t)∈

(0.5,1]×[0,1]×[0,1],

初值条件

u0(x,y)=(x-0.5)sin(x+y),(x,y)∈

[0,0.5]×[0,1],

v0(x,y)=(2x-1)sin(x+y),(x,y)∈

(0.5,1]×[0,1],

边界条件

u(0,y,t)=-0.5e-tsiny,y∈[0,1],

t∈[0,1],

u(x,0,t)=e-t(x-0.5)sin(x),x∈[0,0.5],

t∈[0,1],

u(x,1,t)=e-t(x-0.5)sin(x+1),

x∈[0,0.5],t∈[0,1],

v(1,y,t)=e-tsin(y+1),y∈[0,1],

t∈[0,1],

v(x,0,t)=e-t(2x-1)sin(x),x∈(0.5,1],

t∈[0,1],

v(x,1,t)=e-t(2x-1)sin(x+1),

x∈(0.5,1],t∈[0,1],

界面条件

φ(0.5,y,t)=ψ(0.5,y,t)=0,(y,t)∈

[0,1]×[0,1],

精确解为

u(x,y,t)=e-t(x-0.5)sin(x+y),

(x,y,t)∈[0,0.5]×[0,1]×[0,1],

v(x,y,t)=e-t(2x-1)sin(x+y),

(x,y,t)∈(0.5,1]×[0,1]×[0,1].

构造近似解

U(x,y,t) = (1 -t)((x- 0.5)sin(x+y)) +

(1-x)(-0.5e-tsin(y)-

(1-t)(-0.5sin(y)))+

(1-y)(e-t(x-0.5)sin(x)-

(1-t)((x-0.5)sin(x)))+

y(e-t(x-0.5)sin(x+1)-

(1-t)((x-0.5)sin(x+1))-

(1-x)(-0.5e-tsin(1)-

(1-t)(-0.5sin(1))))+

xy(1-y)tN1(x,y,t,p1),

V(x,y,t)=(1-t)(2x-1)sin(x+y)+

(2x-1)(e-tsin(1+y)-

(1-y)(e-t(2x-1)sin(x)-

(1-t)(2x-1)sin(x)-(2x-1)(e-tsin(1)-

y(e-t(2x-1)sin(x+1)-

(1-t)(2x-1)sin(x+1)-

(2x-1)(e-tsin(2)-(1-t)sin(2))-

(2-2x)(2x-1)y(1-

y)tN2(x,y,t,p2),

其中神经网络N1,N2的隐层神经元数M1=M2=10.

由于本算例空间区域是2维的，若使用Adam 算法求解优化问题(15)，则需选取较多离散训练点，这将导致计算复杂度增加.因此，我们基于随机梯度法思想提出一种针对界面问题的类随机梯度法(算法1) .此算法通过如下步骤更新神经网络N1,N2的权重p1,p2：

1.随机生成一点x=((x0,y0),t0)∈Ω×[0,T],然后在该点计算损失函数关于权重p1,p2的梯度，选择适当步长，沿负梯度方向更新权重;

2.重复步骤1若干次，得到新的权重;

3.将新的权重代入近似解，计算近似解在测试集上的残差ER，判断是否满足停机准则：若满足，则输出优化权重；若不满足，则重复步骤1,2,3.

在算法1中，每次重复步骤1迭代更新权重时，只代入一个随机训练点，对应的离散问题(15)中优化泛函只含一项表达式，其计算复杂度低于多训练点情形.当然，算法1无法避免一般随机算法的共同缺点，即在每次迭代时可能不会按照最优的方向进行梯度下降.但另一方面，由于随机性这类算法存在这样的可能性：跳出当前陷入的局部最优并最终达到全局最优.

算法4.3针对界面问题的一种类随机梯度法 (L-SGD)

输入: 随机生成的初始权重p1,p2,正实数ε,每轮迭代生成总点数I，校正比例β，最大迭代轮数Nmax输出: 最终权重p1,p2,

1: function L-SGD(p1,p2,ε)

2:ER=ε+ 1,n=0

3: whileER>εandn

4:i= 0

5: whilei

6:随机生成一点x=((x0,y0),t0)∈× [0,T]

7:ifx∈Ω1×[0,T] then

8:g=▽p1LossU,确定步长η，p1=p1-ηg

9: else

10: ifx∈Ω2× [0,T] then

11:g=▽p2LossV,确定步长η，p2=p2-ηg

12:else

13:g=▽p1,p2(αLossJ),确定步长η，p=(p1,p2)，p=p-ηg

14: end if

15:end if

16:i=i+1

17:end while

18: - - - 以下j循环是前面i循环的校正，目的是为解决i循环中随机生成的点落在界面上概率较低的问题- - -

19:j=0

20: whilej

21:随机生成界面上一点x=x=((x0,y0),t0) ∈Γ× [0,T]

22:g=▽p1,p2(αLossJ),确定步长η，p=(p1,p2)，p=p-ηg

23:j=j+ 1

24:end while

25: - - - - - - - - - - - - - - -

26:在测试集上计算ER,n=n+1

27:end while

28: return 最终权重p1,p2

29: end function.

在本算例中,优化问题中的参数α=1，神经网络N1,N2的隐层神经元数M1=M2=10，每轮迭代生成总点数I=500，取定步长η=0.001， 校正比例β=6%.

图12 近似解在测试集上的ER误差随迭代轮数变化曲线Fig.12 ERerror curve for the approximate solution computed on the test set

图13 近似解的E2误差随迭代轮数变化曲线

图12和图13分别示出了近似解在测试集上的ER误差和其E2误差随迭代轮数变化的曲线；表2给出不同迭代轮次的近似解在不同时刻的E2误差.图12显示，前10 轮迭代中，ER误差下降较快；迭代20 轮后，其趋于稳定.从图13 和表2可以看出，在前60 轮迭代中，近似解的E2误差下降较快；迭代120 轮后，其趋于稳定.上述结果表明，类随机梯度算法(算法1)是可行的，且能给出较好的数值结果.

表2 近似解不同迭代轮次不同时刻的E2 误差

注5从算例4.1和4.2的数值结果可以看出，当算法迭代到一定次数后，误差趋于稳定.我们分析其原因，主要有以下三个方面：

·算例中神经网络隐层神经元数目较少，导致近似解与真解始终存在一定的误差(参见定理2.1)；

·代入离散优化问题(15)的训练点数有限，不可避免会产生一定的离散误差；

·关于神经网络权重的离散优化问题(15) 是强非线性的，优化算法一般很难得到全局最优解.

本文讨论了求解热传导界面问题的一种单隐层神经网络算法.该算法使用神经网络函数构造满足初边值条件和Dirichlet界面条件的近似解，通过求解由原方程导出的关于神经网络权重的离散优化问题来训练近似解中的单隐层神经网络，使近似解逼近真解.本文也给出了一种基于随机梯度法思想的类随机梯度法来求解相应的离散优化问题. 数值算例验证了神经网络算法的有效性.这一算法既具有近似解表达形式简洁、适用范围广等优点，又易于处理界面条件，便于推广应用于其它类型的界面问题.