1. 引言

科学计算现已随着计算机科学的迅速进展，并肩于理论研究和实验方法，成为科学探索的三大支柱之一[1][2]。在科学计算的众多分支中，求解偏微分方程(PDEs)的数值方法显得尤为关键[3]。许多复杂的物理现象都可通过偏微分方程来建模，然而对这些方程组直接解析求解往往是不可行的[4]。

90年代末，高斯过程(Gaussian process, GP)作为统计学习领域的一种重要工具[5][6]，开始被应用于PDE问题的求解。高斯过程是直接对函数建模的非参数模型，有容易实现、超参数自适应获取，预测值具有概率意义等优点[7]。1992年J.skilling在[8]中提出了用贝叶斯方法求解常微分方程，为后续的研究提供了重要的启示；2011年S. Särkkä在[9]中讨论了高斯过程回归模型的一个扩展，其中测量值被建模为基于高斯过程的线性泛函，估计目标是过程的一般线性算子；2017年M. Raissi在[10]中将偏微分方程描述的物理规律编码为高斯过程先验引入到非参数贝叶斯回归中，在[11]中发展了一种基于参数线性方程的守恒定律表达方法。同年M. Raissi通过使用经典的守恒定律作为结构化先验提高机器学习算法的鲁棒性和效率，提出了使用深度神经网络去计算参数化偏微分方程近似解的PINN方法。2019年P. Rai在[12]中提出了以GP为基函数的PDE参数估计的GPPDE方法，将GP回归模型与状态变量的观测值进行拟合，然后利用GP的导数来获得状态变量的导数，最后调整PDE参数以使GP派生的偏导数满足PDE。该方法适用于扩散方程和Richards方程，且相对Hydrus-1D方法在计算时间上有明显优势；2021年Y. Chen提出TD-GPsolver方法，通过在有限配置点上求解该PDE解，将解近似为高斯过程的最大后验概率(MAP)估计[13]，对椭圆型偏微分方程、Burgers方程和正则化Eikonal方程的实验证明了该方法的有效性。

GPPDE方法和TD-GPsolver方法常被应用在热方程、波动方程等线性偏微分方程中，这些方程由于其数学结构的简洁性使得传统方法可以有效求解。然而，非线性偏微分方程由于解的多样性和非线性动态的存在使得线性方法难以直接适用或效果不佳。本文将GPPDE方法和TD-GPsolver方法扩展到非线性偏微分方程(Allen-Cahn方程和Cahn-Hilliard方程)的求解中，并将两种方法与PINN方法进行比较。

本文的其余部分结构如下：在第2节中，我们详细阐述了基于高斯过程的GPPDE算法和TD-GPsolver算法，并探讨了物理信息神经网络PINN算法；在第3节中我们进行一系列数值实验，将各算法应用到两种非线性偏微分方程中，并对比分析结果，评估了方法的有效性；最后一节对文章进行总结和讨论。

2. 背景知识

2.1. GPPDE方法

假设数据点服从多元高斯分布，其中协方差是其次且平稳的，并且取决于两点之间的距离。对于训练点x和测试点 $x^{\prime }$使用的高斯协方差函数为，

$k\left(x,x^{\prime }\right)={\sigma }_s^2\exp \left(\frac{1}{2}{\displaystyle \sum \frac{{\left(x_j-{x^{\prime }}_j\right)}^2}{l_j^2}}\right)$(1)

其中 $l_j$是第j个自变量的长度，观测值具有独立且同分布的噪声，假设其服从均值为零的正态分布，真实值为 $u_i$，

$\tilde{u}=u_i+N\left(0,{\sigma }_u\right)$(2)

x包含所有自变量 $x_i$， $\tilde{u}$包含所有观测值 ${\tilde{u}}_i$，且服从高斯分布 $u~N\left(m_u\left(x\right),K\right)$， $m_u\left(x\right)$是x的先验均值， $K=k\left(x,x\right)$是协方差矩阵。

使用最大边际似然法[6]用于估计高斯过程的超参数 $\phi =\left[l_s,l_t,{\sigma }_s\right]$，其中 $l_s$， $l_t$分别是空间(x)和时间(t)维度中的长度比例参数，在给定超参数的情况下，边际似然函数通过从真值u和观测值 $\tilde{u}$的联合分布中积分出u获得，

$P\left(\tilde{u}|x,\phi ,{\sigma }_u\right)={\displaystyle \int P\left(\tilde{u},u|x,\phi ,{\sigma }_u\right)\text{d}u}={\displaystyle \int P\left(\tilde{u}|u,x,{\sigma }_u\right)P\left(u|x,\phi \right)\text{d}u}$(3)

使用BFGS算法[14]最大化边际似然函数，完成高斯过程超参数的估计。由于GPPDE方法需要由高斯过程拟合的状态变量关于自变量的导数，导数是一个线性算子，一个高斯过程的导数是另一个高斯过程[15]。对于连续均值的高斯过程一阶导数为，

$p\left({\stackrel{\dot{}}{u}}^{GP}|u,\Phi \right)=p\left({\stackrel{\dot{}}{u}}^{GP}|m_u,\Phi \right)=N\left(\stackrel{\dot{}}{\mu },{\stackrel{\dot{}}{\sigma }}^{GP}\right)$(4)

将偏微分方程以隐函数形式给出，

$f\left(x_1,\cdots ,x_m,u,\frac{\partial u}{\partial x_1},\cdots ,\frac{\partial u}{\partial x_m},\frac{{\partial }^{2}u}{\partial x_1\partial x_1},\cdots ,\frac{{\partial }^{2}u}{\partial x_1\partial x_m},\cdots ,a\right)=0$(5)

其中u是状态变量， $x_j\left(j=1,\cdots ,m\right)$是自变量，a是偏微分方程的一组参数，GPPDE方法对高斯过程的状态值建模，并估计其在偏微分方程中涉及的导数，通过最小化偏微分方程在观测点的残差完成求解。

2.2. TD-GPsolver方法

考虑有界域 $\Omega \subseteq R^d$，其中 $d\ge 1$。非线性偏微分方程给出形式，

$\{\begin{array}{l}P\left(u^{*}\right)\left(x\right)=f\left(x\right),\forall x\in \Omega \\ B\left(u^{*}\right)\left(x\right)=g\left(x\right),\forall x\in \partial \Omega \end{array}$(6)

其中，P是非线性算子，B是边界算子，分别对应f和g。利用在 $\bar{\Omega }$上的有限配置点集合 ${\left\{x_i\right\}}_{i=1}^M$上满足偏微分方程约束的高斯过程逼近真解 $u^{*}$，其中 $x_1,\cdots ,x_{M_{\Omega }}\in \Omega$的内点，边界点为 $x_{M_{\Omega }+1},\cdots ,x_M\in \partial \Omega$，令U为二次Banach空间，相关协方差算子 $K:U^{*}\to U$，考虑最优化问题[16]，

$\{\begin{array}{l}\min \Vert u\Vert \\ \text{s}\text{.t}.P\left(u\right)\left(x_m\right)=f\left(x_m\right),m=1,\cdots ,M\\ B\left(u\right)\left(x_m\right)=g\left(x_m\right),m=M_{\Omega }+1,\cdots ,M\end{array}$(7)

在上述假设下，可以定义泛函 ${\phi }_m^{\left(q\right)}:={\delta }_{x_m}\circ L_q$，同时定义测量向量 $\left[\phi ,u\right]=y\in R^M$也是高斯的，且 $\left[\phi ,u\right]~N\left(0,\Theta \right)$。矩阵 $\Theta$为方阵矩阵，其中 ${\Theta }^{\left(q,j\right)}=k\left({\phi }^{\left(q\right)},{\phi }^{\left(j\right)}\right)$。定义

$y_m=\{\begin{array}{l}f\left(x_m\right),m\in \left\{1,\cdots ,M_{\Omega }\right\}\\ g\left(x_m\right),m\in \left\{M_{\Omega }+1,\cdots ,M\right\}\end{array}$(8)

$F{\left[\phi ,u\right]}_m:=\{\begin{array}{l}P\left({\phi }_m^{\left(1\right)},\cdots ,{\phi }_m^{\left(Q\right)}\right),m\in \left\{1,\cdots ,M_{\Omega }\right\}\\ B\left({\phi }_m^{\left(1\right)},\cdots ,{\phi }_m^{\left(Q_b\right)}\right),m\in \left\{M_{\Omega }+1,\cdots ,M\right\}\end{array}$(9)

求解非线性偏微分方程则转换为考虑优化问题

$\{\begin{array}{l}\min \Vert u\Vert \\ \text{s}\text{.t}.F\left(\left[\phi ,u\right]\right)=y\end{array}$(10)

针对解决无约束优化问题，TD-GPsolver方法采用高斯–牛顿算法的一种有效变体，其中采用自适应块金项与底层核矩阵的离线Cholesky分解一起进行正则化。

2.3. PINN方法

考虑偏微分方程[17]

$\frac{\partial u\left(x,t\right)}{\partial t}+N\left[u\left(x,t\right);\lambda \right]=0$(11)

其中 $u\left(x,t\right)$代表未知函数， $N\left[\cdot ;\lambda \right]$代表由 $\lambda$参数化的非线性微分算子。方程的定义域为 $\Omega$。边界条件初值条件分别是 $B\left(u\left(x,t\right)\right)=0$和 $I\left(u\left(x,0\right)\right)=0$，其中边界条件可以是Dirichlet边界条件，Neumann边界条件，Robin边界条件或混合边界条件。

PINN使用一个被一系列参数化的深度神经网络去近似偏微分方程的未知函数 $u\left(x,t\right)$，即 $u\left(x,t\right)=\tilde{u}\left(x,t\right)$。可以通过最小化方程损失(MSE_f)，边界损失(MSE_b)和初值损失(MSE_i)来学习得到。

3. 数值实验

3.1. 实验设置

为验证前文所述方法的有效性，我们将其应用于两个代表性的非线性偏微分方程：Allen-Cahn方程和Cahn-Hilliard方程。这些方程的参数a均设置为0.01。实验设计包括使用不同规模的离散点数和时间步长，具体配置为：离散点数M= 50与时间步数T= 20；M= 100与T= 40；M= 250与T= 100；以及M= 500与T= 200。在实验设置中，GPPDE方法对观测数据施加了标准差为0.05的噪声，以模拟实际观测条件中可能遇到的数据不确定性。另一方面，TD-GPsolver方法在求解过程中采用了高斯牛顿方法，且迭代步数固定为2次，旨在优化求解过程的效率与精度。

本章节采用平均相对 $L_{2}error$来评估模型的精度，其计算方式如下，

$L_{2}error=\frac{{\displaystyle \sum {\left(u_i-{\stackrel{⌢}{u}}_i\right)}^2}}{N}$

其中 $u_i$代表真实值， ${\hat{u}}_i$的是预测值，N是观测值的数量。

3.2. Allen-Cahn方程

在科学研究领域，Allen-Cahn方程作为一种经典的相场模型，起源于Allen和Cahn关于结晶固体反相位边界运动的描述[18]，已经发展成为解释相变和界面动力学现象的重要非线性偏微分方程。这一方程在材料科学、相变动力学、图像处理以及计算科学等多个学科领域发挥着核心作用，对深入理解物理学中相变过程的本质、界面的演化行为以及相变控制策略提供了理论基础。

考虑Allen-Cahn方程的一维形式如下，

$\frac{\partial u}{\partial t}=a^2\frac{{\partial }^{2}u}{\partial x^2}-u^3+u$(12)

初始条件和边界条件如下，

$\begin{array}{l}u\left(x,0\right)=x^2\cos \left(\pi x\right)\\ u\left(x,t\right)=u\left(-x,t\right)\\ \stackrel{\dot{}}{u}\left(x,t\right)=\stackrel{\dot{}}{u}\left(-x,t\right)\end{array}$(13)

数值实验结果见表1。对于Allen-Cahn方程的计算时间研究，当离散点数为M= 50 (时间步数T= 20)时，GPPDE、TD-GPsolver和PINN方法的计算时间分别为15秒、9秒和600秒。增加离散点至M= 100 (时间步数T= 40)，GPPDE方法的计算时间显著增加至430秒，TD-GPsolver仅需10秒，而PINN的耗时增至670秒。在M= 250 (时间步数T= 100)时，GPPDE方法未能有效求解，TD-GPsolver和PINN的耗时分别为11秒和520秒。当离散点数增至M= 500 (时间步数T= 200)时，GPPDE方法仍未能求解，TD-GPsolver耗时13秒，而PINN耗时700秒。

Table 1.Simulation results of Allem-Cahn equation data

表1. Allem-Cahn方程数据模拟结果

求解精度方面，在离散点数为M= 50 (时间步数T= 20)时，GPPDE、TD-GPsolver和PINN的L2误差分别为0.12、0.333和0.145。将离散点增至M= 100 (时间步数T= 40)，GPPDE、TD-GPsolver和PINN的L2误差分别调整为0.0852、0.0901和0.0638。进一步增加至M= 250 (时间步数T= 100)，GPPDE方法未能求解，TD-GPsolver的L2误差改善至0.0439，PINN的误差最低，为0.00219。在M= 500 (时间步数T= 200)时，GPPDE方法仍未完成求解，TD-GPsolver和PINN的L2误差为0.0235和0.00423。

结果表明针对Allen-Cahn方程，尽管GPPDE在较高离散点数时因计算资源限制无法求解，但其在低离散点设置下展现出可观的求解精度。相较之下，TD-GPsolver方法不仅在低至高的各种离散点数设置中保持了显著的计算效率，而且随着离散点数的增加求解精度显著提高，证明了其在处理大规模和复杂非线性偏微分方程时的强大能力。需要注意的是，PINN方法随着离散点数量增加出现了过拟合现象，这可能是因为数据量与网络容量平衡和优化算法的快速收敛等因素共同作用。

3.3. Cahn-Hilliard方程

Cahn-Hilliard方程，由John W. Cahn和John E. Hilliard于1958年提出[19]，旨在描述两相混合物中自由能密度的演化路径。该方程根植于变分原理，通过对自由能泛函的最小化过程推导出描述体系演化的方程。该方程的核心特征在于其非线性和耗散的双重属性。非线性属性主要源于自由能泛函中的高阶项，这些项负责描绘相界面的复杂曲率效应以及系统的非均质性。而耗散性则通过扩散项实现，负责引导自由能的逐渐降低，促使体系朝向平衡态的方向发展。

$\frac{\partial u}{\partial t}=-\Delta \left(a^2\frac{{\partial }^{2}u}{\partial x^2}-u^3+u\right)$(14)

初始和边界条件如下，

$\begin{array}{l}u\left(x,0\right)=\cos \left(\pi x\right)-\exp \left(-4{\left(\pi x\right)}^2\right)\\ u\left(x,t\right)=u\left(-x,t\right)\\ \stackrel{\dot{}}{u}\left(x,t\right)=\stackrel{\dot{}}{u}\left(-x,t\right)\end{array}$(15)

数值实验结果见表2。对于Cahn-Hilliard方程的计算时间研究，当离散点数为M= 50 (时间步数T= 20)时，GPPDE、TD-GPsolver和PINN方法的计算时间分别为14秒、8秒和160秒。增加离散点至M= 100 (时间步数T= 40)，GPPDE方法的计算时间显著增加至660秒，TD-GPsolver仅需9秒，而PINN的耗时为300秒。在M = 250 (时间步数T= 100)时，GPPDE方法未能有效求解，TD-GPsolver和PINN的耗时分别为11秒和300秒。当离散点数增至M= 500 (时间步数T= 200)时，GPPDE方法仍未能求解，TD-GPsolver耗时12秒，而PINN耗时540秒。

求解精度方面，在离散点数为M= 50 (时间步数T= 20)时，GPPDE、TD-GPsolver和PINN的L2误差分别为0.0214、0.263和0.073。将离散点增至M= 100 (时间步数T= 40)，GPPDE、TD-GPsolver和PINN的L2误差分别调整为0.0117、0.207和0.225。进一步增加至M= 250 (时间步数T= 100)，GPPDE方法未能求解，TD-GPsolver的L2误差改善至0.199，PINN的误差最低，为0.0154。在M= 500 (时间步数T= 200)时，GPPDE方法未能完成求解，TD-GPsolver和PINN的L2误差为0.196和0.0476。

针对Cahn-Hilliard方程，GPPDE方法处理小规模数据集时仍能提供满意的精度，在处理大规模问题上计算资源受限需进行进一步的优化。TD-GPsolver方法凭借其在各种离散点设置下的快速计算能力和较为稳定的精度表现，证明了高斯过程方法是求解非线性偏微分方程的有效工具。PINN方法则在高离散点数设置下展现出较优的求解精度，但对计算资源和计算时间要求较高。

Table 2.Simulation results of Cahn-Hilliard equation data

表2.Cahn-Hilliard方程数据模拟结果

4. 结论与展望

本文重点分析了两种基于高斯过程的数值解法(TD-GPsolver方法和GPPDE方法)在处理复杂NPDEs中的应用。通过对Allen-Cahn方程和Cahn-Hilliard方程的求解进行数值实验，我们评估了这两种方法在不同离散化水平下的求解精度(以L2范数误差衡量)和计算效率(以计算时间为指标)，并以PINN方法作为对比。TD-GPsolver方法利用时间动态的高斯过程优化了计算过程，并在维持低误差率的同时显著提高了求解速度，表现出在处理大规模数据集时拥有的高效率和稳定性。另一方面，GPPDE虽然在高离散点数量下表现不佳，但在处理小规模数据集时仍能提供较高精度，特别是在对概率分布有详细要求的应用场景中。PINN方法需要大量计算资源和计算时间，适合应用精度要求极高的场景中。

总结来说，TD-GPsolver和GPPDE方法都体现了高斯过程在数值求解非线性偏微分方程中的独特价值和潜力。未来研究可以进一步探索TD-GPsolver方法在其他类型的偏微分方程求解中的应用，以及如何通过算法优化进一步提升其在低离散化水平下的精度表现。此外，结合TD-GPsolver的高效性与PINN的高精度特性，开发新的算法以应对更广泛的科学和工程问题，也是未来研究的重要趋势。

致 谢

作者感谢主编、编委会和匿名审稿人提出许多宝贵意见，这些意见对促进文章有重大改进。作者对撰写文章时提供帮助的师生表示感谢。

参考文献