求解时间分数阶相场微分方程的自适应分数阶物理信息网络

Adaptive Fractional Physical Information Neural Network for Solving Time Fractional Phase Field Partial Differential Equations

摘要: 本文提出具有自适应权重的分数阶物理信息神经网络(adaptive-fPINN-PQI)求解时间分数阶偏微分方程。首先，利用Hadamard有限部分积分意义上的分段二次插值(PQI)对时间分数阶导数进行离散。其次，为降低自动微分引入的误差，本文采用中心差分法代替自动微分求导，计算空间各阶偏导数，提高了预测解精度。此外，本文构建的自适应权重残差网络，基于残差网络架构有效防止梯度消失。并通过建立自适应权重来自动调整不同损失项的权重，显著平衡其梯度，进一步提升预测解精度。最后，将adaptive-fPINN-PQI用于求解时间分数阶相场偏微分方程，证明了该网络的高精度和高效率。

Abstract: In this paper, an accurate fractional physical information neural network with an adaptive learning rate (adaptive-fPINN-PQI) is proposed for solving time fractional partial differential equations. Firstly, the time-fractional derivative in the Caputo sense is discretized by piecewise quadratic interpolation (PQI) in the sense of the Hadamard finite-part integral. Next, the central difference method is used instead of automatic differentiation to reduce the error of automatic differentiation. Besides, the present adaptive-fPINN-PQI is based on the ResNet architecture to effectively overcome the issue of gradient vanishing. The adaptive learning rate is constructed to automatically adjust the weights of different loss terms, significantly balancing their gradients and improving the accuracy of the predicted solutions. Finally, time fractional phase field equations have been solved using the proposed adaptive-fPINN-PQI to demonstrate its high precision and efficiency.

1. 引言

分数阶微积分的非局域性和记忆性使其更能描述复杂的数学物理问题，因此许多学者对该问题进行了大量的研究。目前，已经发展出一系列有效的求解分数阶微分方程的数值方法，如有限差分法 [1] 、有限元法 [2] 、谱方法 [3] 、无网格法(MFM) [4] 等。其中，T.A.M.和B.I. [5] 提出了一种一阶有限差分(L1)的时间离散方案。Deng [6] 研究了分数阶Fokker-Planck方程的有限元方法。Li等 [7] 提出了一种无网格数值逼近一维问题的方法——“有限积分法(FIM)”。Li等 [8] 证明了FIM在求解多维问题时具有较高的稳定性、效率和精度。

深度学习作为一种求解整数阶偏微分方程的数值技术已广泛应用。许多基于数据驱动的深度神经网络被提出求解偏微分方程 [9] 。Wang等 [10] 将残差神经网络与数据驱动方法相结合，有效地解决了高维问题。Yohai等 [11] 提出了一种数据驱动的离散化方案，利用神经网络估计空间导数，并对空间导数进行端到端优化。然而，数据驱动的网络无法仅通过物理约束来解决这个问题，因为它们需要大量的标记数据。针对数据驱动的局限性，研究人员提出物理信息神经网络 [12] 。PINNs (物理信息神经网络)将偏微分方程作为惩罚项加入到损失函数中，使网络能够学习物理方程中隐藏的信息。目前，求解经典偏微分方程的物理信息神经网络技术得到了很大的发展，并在实际应用中显示出良好的效果 [13] [14] [15] 。然而，对于整数阶偏微分方程，利用网络的自动微分可以很容易地推导出空间和时间的任意阶导数，而分数阶导数无法通过自动微分获得，因此构建稳定高效的分数阶偏微分方程神经网络已成为深度学习领域的研究热点。fPINN [16] 使用L1方案离散时间分数阶导数，利用GL方法离散空间拉普拉斯分数阶算子。Rostami [17] 将级数展开作为替代解，实现了高阶线性分数阶偏微分方程(FPDEs)的更精确解析。Wang等 [18] 将分数阶物理信息与神经网络的非线性函数相结合，提高了时间分数阶相场模型的建模和求解效果。

经典的一阶有限差分(L1)格式具有良好的稳定性和数值精度。但由于其收敛阶较低，可能需要更多的网格节点才能获得与高阶差分格式相同的数值精度。Li等 [19] 提出的Hadamard有限部分积分意义上的分段正交插值(PQI)具有三阶收敛性，可以提高预测解的精度，提高模型的效率。

因此，本文在时间离散化中采用PQI格式，快速减小控制方程的残差损失，得到更精确的解。为了平衡损失函数中不同损失项的梯度，在损失函数中提出自适应权重，防止网络陷入局部最优解，大大提高了网络的可训练性、收敛性和精度。在此基础上，采用中心差分法代替网络中的自动微分来计算空间各阶偏导数以进一步提高系统精度。

2. 数学模型

2.1. 问题描述

本文考虑时间分数阶偏微分方程的如下形式：

$\begin{array}{l}{}_0{}^{C}D{}_t^{\alpha }u\left(\mathbb{X},t\right)-Lu\left(\mathbb{X},t\right)=f\left(\mathbb{X},t\right),\mathbb{X}\in \Omega ,t\in \left[0,T\right],\\ u\left(\mathbb{X},0\right)=g_1\left(\mathbb{X}\right),\mathbb{X}\in \Omega ,\\ u\left(\mathbb{X},t\right)=g_2\left(\mathbb{X}\right),\mathbb{X}\in \partial \Omega ,t\in \left[0,T\right],\end{array}$ (1)

其中， ${}_0{}^{C}D{}_t^{\alpha }u\left(\mathbb{X},t\right)$ 是 $u\left(\mathbb{X},t\right)$ 对时间t的卡普托意义分数阶导数， $\alpha$ 是导数的阶数，L表示微分算子， $\Omega$ 是问题域， $\partial \Omega$ 是 $\Omega$ 的边界。这里g₁和g₂是两个已知的连续函数。

2.2. 时间分数阶导数的离散格式

分数阶导数在卡普托意义上定义为

${}_0{}^{c}D{}_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-x\right)}^{-\alpha }{f}^{\prime }\left(t\right)\text{d}x},0<\alpha <1,t>0,$ (2)

其中Γ表示伽马函数， $\alpha$ 是时间分数阶导数的阶数， $\alpha \in \left(0,1\right)$ 。区间[0,T]均匀划分为 $0=t_0<t_1<t_2<\cdots <t_{2s}<t_{2s+1}<\cdots <t_{2N_T}=T$ ，其中s和N_T表示位置整数。

在分段求积分形式中，首先使用黎曼-刘维尔分数阶导数表示式(2)为

$D_t^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(\alpha \right)}\hat{\int }{\left(t-x\right)}^{-1-\alpha }f\left(t\right)\text{d}x.$ (3)

基于分段二次插值多项式，有：

$D_t^{\alpha }f\left(\mathbb{X},t_{2s}\right)=\Delta t^{-\alpha }{\displaystyle \underset{k=0}{\overset{2s}{\sum }}{\nu }_{k,2s}\hat{f}\left(\mathbb{X},t_{2s-k}\right)}+R_2^{2s},$ (4)

式中 $\left|R_2^{2s}\right|\le C\Delta t^{3-\alpha }\underset{0<t<1}{\max }\left|{\stackrel{\dot{}}{\hat{f}}}^{\left(3\right)}\left(\mathbb{X},t_{2s}\right)\right|$ ， ${\stackrel{\dot{}}{\hat{f}}}^{\left(3\right)}$ 为三阶导数，C为有限常数。权重 ${\nu }_{k,2s}$ ， $k=0,1,2,\cdots ,2s$ ， $s=1,2,\cdots ,N_T$ 具有如下形式：

${\nu }_{k,2s}=\frac{{\mathbb{W}}_{k,2s}}{\Gamma \left(3-\alpha \right)},k=1,2,\cdots ,2s,$ (5)

其中

${\mathbb{W}}_{k,2s}=\{\begin{array}{l}{2}^{-\alpha }\left(2-\alpha \right),\text{for}k=0,\hfill \\ 2^{2-\alpha }\left(-\alpha \right),\text{for}k=1,\hfill \\ \frac{1}{2}{W}_0\left(2\right)+\left(2^{-\alpha }\right){\alpha }^2,\text{for}k=2,\hfill \\ -W_1\left(m\right),\text{for}k=2m-1,m=2,3,\cdots ,s,\hfill \\ \frac{1}{2}\left(W_2\left(m\right)+W_0\left(m+1\right)\right),\text{for}k=2m,m=2,3,\cdots ,s-1,\hfill \\ \frac{1}{2}{W}_2\left(s\right),\text{for}k=2s.\hfill \end{array}$ (6)

$W_i\left(k\right)$ 和 $V_i\left(k,\alpha \right)$ 的计算公式为：

$\begin{array}{l}{W}_i\left(k\right)=V_i\left(k,\alpha \right)\left({\left(2k\right)}^{-\alpha }-{\left(2k-2\right)}^{-\alpha }\right)\left(1-\alpha \right)\left(2-\alpha \right)-Z_i\left(k,\alpha \right)\left(-\alpha \right)\left(2-\alpha \right)\\ +\left({\left(2k\right)}^{2-\alpha }-{\left(2k-2\right)}^{2-\alpha }\right)\left(-\alpha \right)\left(1-\alpha \right),i=0,1,2\\ V_i\left(k,\alpha \right)=\{\begin{array}{l}\left(2k-1\right)\left(2k\right),i=0,\hfill \\ \left(2k-2\right)\left(2k\right),i=1,\hfill \\ \left(2k-1\right)\left(2k-2\right),i=2,\hfill \end{array}\\ Z_i\left(k,\alpha \right)=\{\begin{array}{l}\left(4k-1\right)\left({\left(2k\right)}^{1-\alpha }-{\left(2k-2\right)}^{1-\alpha }\right),i=0,\hfill \\ \left(4k-2\right)\left({\left(2k\right)}^{1-\alpha }-{\left(2k-4\right)}^{1-\alpha }\right),i=1,\hfill \\ \left(4k-3\right)\left({\left(2k\right)}^{1-\alpha }-{\left(2k-4\right)}^{1-\alpha }\right),i=2.\hfill \end{array}\end{array}$ (7)

这样可以将黎曼–刘维尔分数阶导数转化为卡普托意义上的分数阶导数，形式如下：

${}_0{}^{c}D{}_t^{\alpha }\hat{f}\left(\mathbb{X},t_{2s}\right)=D_t^{\alpha }\left(\hat{f}\left(\mathbb{X},t_{2s}\right)-\hat{f}\left(\mathbb{X},t_0\right)\right)=\Delta t^{-\alpha }{\displaystyle \underset{k=0}{\overset{2s}{\sum }}{\tilde{\nu }}_{k,2s}}\hat{f}\left(\mathbb{X},t_{2s-k}\right)+R_2^{2s},$ (8)

其中， ${\tilde{\nu }}_{k,2s}=\{\begin{array}{l}{\nu }_{k,2s},k=0,1,2,\cdots ,2s-1\hfill \\ {\nu }_{k,2s}-\frac{{\left(2s\right)}^{-\alpha }}{\Gamma \left(1-\alpha \right)},k=2s\hfill \end{array}$ 。

2.3. 空间偏导数的离散格式

为了更精确地近似偏导数，避免网络自动微分引入的误差，本文采用中心差分来近似空间的各阶偏导数，其计算公式为：

在空间近似一阶导数 $\frac{\partial u\left(x\right)}{\partial x}$ 时，本文采用的差分格式为：

$\frac{\partial u\left(x\right)}{\partial x}=\frac{{\hat{u}}_e-{\hat{u}}_w}{\Delta x},$ (9)

${\hat{u}}_e$ 和 ${\hat{u}}_w$ 分别为( $x+\frac{\Delta x}{2}$ )和( $x-\frac{\Delta x}{2}$ )位置处的导数， $\Delta x$ 应为常规数值格式中相邻两点之间的距离。相应的导数项可推导为：

$\frac{\partial u\left(x\right)}{\partial x}\cong {\frac{\partial u\left(x\right)}{\partial x}|}_{cd2}=\frac{{{\hat{u}}_e|}_{cd2}-{{\hat{u}}_w|}_{cd2}}{\Delta x}=\frac{\hat{u}\left(x+\Delta x;w\right)-\hat{u}\left(x-\Delta x;w\right)}{2\Delta x}$ (10)

3. 网络架构

fPINNs (分数阶物理信息神经网络) [16] 利用一阶有限差分(L1)格式来近似时间分数阶导数，从而使分数阶偏微分方程的残差 ${\displaystyle \underset{k=0}{\overset{2s}{\sum }}{\varphi }_{k,2s}u\left(x,y,t;\theta \right)}+L\tilde{u}\left(x,y,t;\theta \right)-f\left(x,y,t\right)$ 可以作为损失项加入到损失函数中，使网络能够更准确学习物理方程中隐藏的信息，其网络结构如图1所示。

Figure 1. The schematic of fPINN

图1. fPINN网络结构

fPINN中不同损失项的系数通常都固定为1，如图1所示。然而，由于不同损失项的梯度不平衡，会导致网络过于关注边界条件、初始条件或控制方程，从而陷入局部最优。为了解决损失项梯度不平衡问题，本文提出自适应权重来调整不同损失项的梯度，保证网络在反向传播过程中关注到所有损失项。

本文的损失函数定义为残差损失项、边界损失项和初始损失项的均方误差(MSE)：

$\begin{array}{l}L\left(\theta \right)=L_r\left(\theta \right)+{\lambda }_i{L}_{u_i}\left(\theta \right)+{\lambda }_b{L}_{u_b}\left(\theta \right),\\ L_r\left(\theta \right)=\frac{1}{N_f}{\displaystyle \underset{i=1}{\overset{N_f}{\sum }}{\left|{\displaystyle \underset{k=0}{\overset{2s}{\sum }}{\varphi }_{k,2s}\tilde{u}\left({\mathbb{X}}_i,t_{2s-k}\right)-f_i}\right|}^2},\\ L_{u_i}\left(\theta \right)=\frac{1}{N_0}{\displaystyle \underset{i=1}{\overset{N_0}{\sum }}{\left|\tilde{u}\left({\mathbb{X}}_i,0\right)-g_{1_i}\right|}^2},\\ L_{u_b}\left(\theta \right)=\frac{1}{N_b}{\displaystyle \underset{i=1}{\overset{N_b}{\sum }}{\left|\tilde{u}\left({\mathbb{X}}_i,t_i\right)-g_{2_i}\right|}^2}.\end{array}$ (11)

其中 ${\lambda }_b$ ， ${\lambda }_i$ 分别表示边界损失项和初始损失项的自适应权重，θ为包含权重和偏置的网络参数，使用梯度下降法更新和迭代：

${\theta }_{n+1}={\theta }_n-\eta {\nabla }_{\theta }L\left(\theta \right)={\theta }_n-\eta \left({\nabla }_{\theta }{L}_r\left(\theta \right)+{\nabla }_{\theta }{L}_{u_i}\left(\theta \right)+{\nabla }_{\theta }{L}_{u_b}\left(\theta \right)\right),$ (12)

式中η为学习率， ${\nabla }_{\theta }{L}_r\left(\theta \right),{\nabla }_{\theta }{L}_{u_i}\left(\theta \right),{\nabla }_{\theta }{L}_{u_b}\left(\theta \right)$ 分别为残差损失项、边界损失项和初始损失项的梯度值。自适应权重需要满足：

$\begin{array}{l}\forall \epsilon >0,\left|{\nabla }_{\theta }{L}_r\left(\theta \right)-{\lambda }_{u_i}{\nabla }_{\theta }{L}_{u_i}\left(\theta \right)-{\lambda }_{u_b}{\nabla }_{\theta }{L}_{u_b}\left(\theta \right)\right|<\epsilon ,\\ 其中{\lambda }_{u_i}=\frac{\overline{\left|{\nabla }_{\theta }{L}_r\left(\theta \right)\right|}}{\overline{\left|{\nabla }_{\theta }{L}_{u_i}\left(\theta \right)\right|}},{\lambda }_{u_b}=\frac{\overline{\left|{\nabla }_{\theta }{L}_r\left(\theta \right)\right|}}{\overline{\left|{\nabla }_{\theta }{L}_{u_b}\left(\theta \right)\right|}}\end{array}$ (13)

其中 $\overline{\left|{\nabla }_{\theta }{L}_r\left(\theta \right)\right|},\overline{\left|{\nabla }_{\theta }{L}_{u_i}\left(\theta \right)\right|}$ 和 $\overline{\left|{\nabla }_{\theta }{L}_{u_b}\left(\theta \right)\right|}$ 分别表示残差损失项、边界损失项和初始损失项相对于参数θ的梯度均值。通过大量的数值实验，作者发现当 ${\lambda }_{u_i}=\frac{\max \left|{\nabla }_{\theta }{L}_r\left(\theta \right)\right|}{\overline{\left|{\nabla }_{\theta }{L}_{u_i}\left(\theta \right)\right|}}$ 和 ${\lambda }_{u_b}=\frac{\max \left|{\nabla }_{\theta }{L}_r\left(\theta \right)\right|}{\overline{\left|{\nabla }_{\theta }{L}_{u_b}\left(\theta \right)\right|}}$ 时，更有利于网络的损失收敛。

此外，考虑到ResNet [20] (图2)可以有效解决fPINN中使用的传统全连接神经网络(DNN)的梯度消失问题，因此，在本文中，从第二层隐藏层开始采用ResNet。

Figure 2. The architecture of ResNet block

图2. 残差块结构

用公式(8) (10)离散式(1)，结合边界条件残差和初始条件残差构成物理信息网络的损失函数。并将本文提出的自适应权重和现有的ResNet块应用到网络结构中，得到的自适应fPINN网络结构如图3所示。

Figure 3. The adaptive-fPINN architecture

图3. 自适应fPINN网络结构

4. 数值实验与结果

4.1. 对比网络

在本节中，将本文提出的网络应用到求解1D时间分数阶相场偏微分方程和2D时间分数阶相场偏微分方程。为了更好地说明网络的精度和性能，将与fPINN [20] 结果进行对比。在网络参数方面，本文提出的网络(adaptive-fPINN)与对比网络采用完全相同的参数设置，具体为：使用Adam算法对损失函数进行优化，Xavier用于初始化，神经网络的学习率为1 × 10^-3，隐藏层数设置为5，神经元数和迭代次数分别固定为200和6 × 10⁴，激活函数使用tanh。

4.2. 1D时间分数阶相场方程

本文对具有Dirichlet边界条件的1D时间分数阶相场偏微分方程进行了研究，方程定义为：

$\begin{array}{l}{D}_t^{\alpha }u\left(\mathbb{X},t\right)-{\epsilon }^2\Delta u-u+u^3=f\left(\mathbb{X},t\right),\mathbb{X}\in \Omega ,t\in \left[0,T\right],\\ u\left(\mathbb{X},0\right)=g_1\left(\mathbb{X}\right),\mathbb{X}\in \Omega ,\\ u\left(\mathbb{X},t\right)=g_2\left(\mathbb{X}\right),\mathbb{X}\in \partial \Omega ,t\in \left[0,T\right].\end{array}$ (14)

其中 $\Delta$ 表示为拉普拉斯算子， $\Omega ={\left[0,1\right]}^d$ ， $d=1,2$ 并且 $\partial \Omega$ 是 $\Omega$ 的边界，固定 $\epsilon =1,T=1$ 。

考虑满足该方程的一个精确解 $u\left(x,t\right)=\left(1+t^2\right)\sin \left(2\pi x\right)$ ，根据式(14)，可得：

$f\left(x,t\right)=\frac{\Gamma \left(3\right)}{\Gamma \left(3-\alpha \right)}{t}^{2-\alpha }\sin \left(2\pi x\right)+4{\pi }^2\left(1+t^2\right)\sin \left(2\pi x\right)-\left(1+t^2\right)\sin \left(2\pi x\right)+{\left[\left(1+t^2\right)\sin \left(2\pi x\right)\right]}^3$ (15)

分别利用fPINN和adaptive-fPINN求解该问题，结果见表1。通过比较表中的数据，可明显观察到本文所提出的adaptive-fPINN的误差比fPINN的误差小1~2个数量级，说明其可显著地提高网络求解精度。此外，图4绘制了这两种模型在迭代过程中，预测解的相对误差。从该箱体图中可以看出，adaptive-fPINN模型在整个迭代过程中具有更短的箱长，表明所提出的网络更稳定。而且adaptive-fPINN的误差始终小于fPINN的误差，总体说明adaptive-fPINN的解更精确、更稳定。

Table 1. Comparison of the relative errors of the prediction solutions using different networks for the 1D time fractional AC equation when t = 1

表1. t = 1时一维时间分数相场方程网络预测解的相对误差比较

4.3. 2D时间分数阶相场方程

在这里考虑具有精确解为 $u\left(x,t\right)=\left(1+t^2\right)\sin \left(2\pi x\right)\sin \left(2\pi y\right)$ 的2D时间分数阶相场偏微分方程，根据式(14)，可得源项为：

$\begin{array}{l}f\left(x,y,t\right)=\frac{\Gamma \left(3\right)}{\Gamma \left(3-\alpha \right)}{t}^{2-\alpha }\sin \left(2\pi x\right)\sin \left(2\pi y\right)+8{\pi }^2\left(1+t^2\right)\sin \left(2\pi x\right)\sin \left(2\pi y\right)\\ -\left(1+t^2\right)\sin \left(2\pi x\right)\sin \left(2\pi y\right)+{\left[\left(1+t^2\right)\sin \left(2\pi x\right)\sin \left(2\pi y\right)\right]}^3\end{array}$ (16)

Figure 4. The relative errors of the solutions against the number of iterations by various models at $\alpha =0.8$

图4. $\alpha =0.8$ 时，网络迭代过程中不同模型解的相对误差

神经网络fPINN和adaptive-fPINN求解该问题的预测解相对误差如表2所示。结果表明，adaptive-fPINN-PQI模型解的精度比fPINN平均高一个数量级。

Table 2. Comparison of the relative errors of the prediction solutions using different networks for the 2D time fractional AC equation when t = 1

表2. t = 1时二维时间分数相场方程网络预测解的相对误差比较

图5绘制了fPINN和adaptive-fPINN解及其误差分布。从图中可以看出，本文提出的网络预测性能要优于对比网络，在整个预测区域内预测精度降低大约1个数量级，由此可见所提出的自适应fPINN模型在求解时间分数阶方程方面表现优异。

(a) fPINN(b) adaptive-fPINN

Figure 5. The distribution of predicted solutions and errors of fPINN and adaptive-fPINN

图5. fPINN和adaptive-fPINN网络的预测解分布和绝对误差分布

5. 结论

本文构造了一个求解时间分数阶偏微分方程的自适应fPINN网络。该结构的自适应权重通过有效平衡不同损失项的梯度，使所有损失项都能快速下降。所提出的神经网络采用了分段二次插值的高阶格式，使得网络损失可以快速下降，提高了训练速度和效率。在此基础上利用中心差分法代替物理信息网络中的自动微分技术，一定程度上降低了近似误差。最后通过求解两个时间分数阶相场偏微分方程实验有效验证了所提出模型的精度、收敛性和稳定性。在今后的工作中，将进一步研究不同定义下的分数阶导数，特别是解决Hadamard型分数阶问题，开发一个更广泛适用的机器学习模型解决分数阶问题。

基金项目

山西省青年科学研究基金项目(No. 20210302124159)；国家自然科学基金项目(No. 12102283)。

NOTES

^*通讯作者。