【PINN】DeepXDE学习训练营(16)——operator-diff_rec_aligned

本文链接：https://blog.csdn.net/2303_77200324/article/details/147677494

一、引言

随着人工智能技术的飞速发展，深度学习在图像识别、自然语言处理等领域的应用屡见不鲜，但在科学计算、工程模拟以及物理建模方面，传统的数值方法仍然占据主导地位。偏微分方程（Partial Differential Equations, PDEs）作为描述自然界中众多复杂现象的重要数学工具，在物理、化学、工程、金融等领域具有广泛应用。然而，伴随着高维度、多变量、复杂边界条件等挑战，传统数值求解方法面临效率低、适应性差等困境。

近年来，深度学习的崛起为科学计算带来了全新的解决思路。其中，以深度偏微分方程（Deep PDE）为代表的研究方向，通过结合神经网络与偏微分方程的理论，成功开发出高效、灵活的求解方案。这种方法不仅可以克服传统方法的局限，还能应对高维、复杂几何等问题。

作为深度偏微分方程领域的开源工具库，DeepXDE（Deep Learning for Differential Equations）由lululxvi团队精心开发，凭借其强大的功能、易用的接口和丰富的示例，受到学术界与工业界的广泛关注。本文将系统介绍DeepXDE的基本内容与应用价值，深入探讨其核心技术原理，分享环境配置与运行技巧，并结合实际案例进行分析，最后对未来发展趋势进行展望。

二、DeepXDE的用途

DeepXDE，一个基于TensorFlow和PyTorch的深度学习微分方程求解库，应运而生。它提供了一个简洁、高效且易于使用的框架，使得研究人员和工程师能够利用深度学习技术求解各种类型的微分方程，包括常微分方程（ODEs）、偏微分方程（PDEs）、积分微分方程（IDEs）以及分数阶微分方程（FDEs）。

DeepXDE旨在提供一站式的深度学习框架，用于高效求解各种偏微分方程，包括但不限于：

1. 传统偏微分方程求解

定常和非定常问题：热传导方程、波动方程、拉普拉斯方程、扩散方程等。
线性和非线性方程：支持线性边界条件，也能处理非线性、非局部问题。

2. 高维偏微分方程

在高维空间中，传统数值方法面临“维数灾难”。DeepXDE利用神经网络天然的高维逼近能力，有效解决高维PDE，如贝尔曼方程、多体问题等。

3. 复杂几何和边界条件

支持任意复杂的几何区域、非均匀边界条件，极大扩展了求解的适用范围。

4.参数逆问题和数据驱动建模

整合数据，使模型在已知部分信息的基础上进行参数识别、反演问题求解。

5. 动态系统和时间演化

支持带有时间变量的演化问题，模拟动态过程。

6. 结合有限元、有限差分等方法

虽然核心为神经网络，但兼容各种数值方法，提供灵活的求解策略。

7. 教育科研与工程实践

丰富的案例与接口帮助科研人员快速验证理论，工程师实现快速设计优化。

总结而言，DeepXDE不仅是一个纯粹的数学工具，更是工程实践中的“聪明助手”，帮助用户以信赖深度学习的方式突破传统技术瓶颈，实现创新性的科学计算。

三、核心技术原理

DeepXDE的核心思想是利用神经网络作为逼近器，通过构造损失函数，使网络能在满足偏微分方程边界条件的前提下逼近真实解。以下详细阐释其原理基础。

1. 神经网络逼近偏微分方程解

假设待求解的偏微分方程可以写成：

$\mathcal{N}[u](x)=0,x\in\Omega,$

配合边界条件

$\mathcal{B}[u](x)=g(x),x\in\partial \Omega,$

这里， $\mathcal{N}$ 代表微分算子， $\mathcal{B}$ 代表边界条件算子。

DeepXDE利用深度神经网络 𝑢𝜃(𝑥) 作为解的逼近，参数为 𝜃 。通过自动微分（AutoDiff），网络可以自然求出 𝑢𝜃 的各阶导数，从而在网络定义的每个点上计算微分方程的残差。

2. 损失函数设计

训练模型的目标是最小化残差，使神经网络逼近满足偏微分方程的解。损失函数由两部分组成：

方程残差部分：

$\mathcal{L}_{\mathit{PDE}}=\frac{1}{N_{f}}\sum_{i=1}^{N_{f}}|\mathcal{N}[u_{ \theta}](x_{f}^{(i)})|^{2},$

其中， $x{_{f}}^{(i)}$ 为采样点，用于评估微分残差。

边界条件部分：

$\mathcal{L}_{\mathit{BC}}=\frac{1}{N_{b}}\sum_{i=1}^{N_{b}}|\mathcal{B}[u_{ \theta}](x_{b}^{(i)})-g(x_{b}^{(i)})|^{2},$

结合整体目标函数：

$\mathcal{L}(\theta)=\lambda_{f}\mathcal{L}_{\mathit{PDE}}+\lambda_{b} \mathcal{L}_{\mathit{BC}},$

这里 $\lambda {_{f}}$ 、 $\lambda {_{b}}$ 为调节系数。

3. 自动微分（AutoDiff）技术

深度学习框架如TensorFlow或PyTorch提供自动微分功能，方便快速计算神经网络输入的微分，自动应用链式法则求导，极大简化偏微分方程的数值差分表达。

4. 训练优化方法

利用成熟的梯度下降（SGD）、Adam等优化算法，通过反向传播调节神经网络参数，使损失函数达到最小。

5. 样本生成和采样策略

采样点生成：采用随机采样、拉丁超立方（Latin Hypercube Sampling）或网格采样来选取训练点。
自适应采样：在训练过程中，根据误差分布调整采样点，提高训练效率。

6. 复杂边界与几何的处理

采用非结构化的几何描述和SDF（Signed Distance Function）结合，保证不同几何形状的灵活支持。

7. 逆问题与数据融合

在已知数据集上引入数据损失，使模型不仅满足PDE，也通过端到端训练实现数据匹配，增强实际适用性。

五、代码详解

"""支持的后端：tensorflow.compat.v1、tensorflow、pytorch、paddle"""
import deepxde as dde
import matplotlib.pyplot as plt
import numpy as np

from ADR_solver import solve_ADR  # 导入自定义的求解随机偏微分方程的函数

# 定义偏微分方程（PDE）
def pde(x, y, v):
    D = 0.01  # 扩散系数
    k = 0.01  # 非线性系数
    dy_t = dde.grad.jacobian(y, x, j=1)  # 计算y关于t的一阶偏导
    dy_xx = dde.grad.hessian(y, x, j=0)   # 计算y关于x的二阶偏导（拉普拉斯算子）
    # PDE：时间导数 - 扩散项 + 非线性项 - 控制变量v
    return dy_t - D * dy_xx + k * y**2 - v

# 定义空间区间
geom = dde.geometry.Interval(0, 1)
# 定义时间区间
timedomain = dde.geometry.TimeDomain(0, 1)
# 建立空间-时间几何区域的组合
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

# 设置边界条件（边界上的值为0）
bc = dde.icbc.DirichletBC(geomtime, lambda _: 0, lambda _, on_boundary: on_boundary)
# 设置初值条件（初始时刻的值为0）
ic = dde.icbc.IC(geomtime, lambda _: 0, lambda _, on_initial: on_initial)

# 构建时间偏微分方程（TimePDE）对象
pde = dde.data.TimePDE(
    geomtime,
    pde,
    [bc, ic],
    num_domain=200,    # 训练中采样的空间和时间点数（内部）
    num_boundary=40,   # 边界点数
    num_initial=20,     # 初始条件点数
    num_test=500,       # 测试点数
)

# 定义特征空间（高斯随机场）
func_space = dde.data.GRF(length_scale=0.2)

# 评估点：在x方向采样50个点
eval_pts = np.linspace(0, 1, num=50)[:, None]

# 构建PDE操作数据，支持操作（PDE和特征空间的乘积）
# 这里function_variables=[0]代表只对第一个变量进行操作
data = dde.data.PDEOperatorCartesianProd(
    pde, func_space, eval_pts, 1000, function_variables=[0], num_test=100, batch_size=50
)

# 构建DeepONet模型（支持操作的多输入结构）
net = dde.nn.DeepONetCartesianProd(
    [50, 128, 128, 128],  # 输入特征路径网络结构
    [2, 128, 128, 128],   # 树路径网络结构（输出空间的维度为2（x,t））
    "tanh",               # 激活函数
    "Glorot normal",      # 权重初始化方式
)

# 创建模型
model = dde.Model(data, net)
# 编译模型，优化器为Adam，学习率0.0005
model.compile("adam", lr=0.0005)
# 训练模型（训练20,000次）
losshistory, train_state = model.train(iterations=20000)
# 绘制训练损失变化
dde.utils.plot_loss_history(losshistory)

# 生成特征空间的采样特征（随机样本）
func_feats = func_space.random(1)
# 生成x轴样本点（0到1之间）
xs = np.linspace(0, 1, num=100)[:, None]
# 在这些点上计算对应的随机控制函数值
v = func_space.eval_batch(func_feats, xs)[0]

# 调用自定义的求解器，获得真实的u(x,t)
x, t, u_true = solve_ADR(
    0, 1, 0, 1,                    # 空间和时间区间
    lambda x: 0.01 * np.ones_like(x),  # 扩散系数
    lambda x: np.zeros_like(x),        # 初始条件（非变化）
    lambda u: 0.01 * u**2,             # 非线性项
    lambda u: 0.02 * u,                  # 非线性系数
    lambda x, t: np.tile(v[:, None], (1, len(t))),  # 控制变量（随机控制函数）
    lambda x: np.zeros_like(x),            # 边界条件
    100,                                 # 空间点数
    100,                                 # 时间点数
)

# 真实解（二维数组，转置，使得行对应时间，列对应空间）
u_true = u_true.T

# 绘制真实解的热力图
plt.figure()
plt.imshow(u_true)
plt.colorbar()

# 生成控制变量的批量（随机特征）
v_branch = func_space.eval_batch(func_feats, np.linspace(0, 1, num=50)[:, None])

# 创建空间-时间网格，用于绘制预测结果
xv, tv = np.meshgrid(x, t)
# 样本点（x,t）扁平化为二维数组
x_trunk = np.vstack((np.ravel(xv), np.ravel(tv))).T
# 用训练好的模型进行预测
u_pred = model.predict((v_branch, x_trunk))
# 预测结果reshape为（空间点数，时间点数）
u_pred = u_pred.reshape((100, 100))

# 计算相对L2误差
print(dde.metrics.l2_relative_error(u_true, u_pred))
# 绘制