【PINN】DeepXDE学习训练营(34)——pinn_forward-fractional_Poisson

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

一、引言

随着人工智能技术的飞速发展，深度学习在图像识别、自然语言处理等领域的应用屡见不鲜，但在科学计算、工程模拟以及物理建模方面，传统的数值方法仍然占据主导地位。偏微分方程（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，也通过端到端训练实现数据匹配，增强实际适用性。

五、代码详解

"""Backend supported: tensorflow.compat.v1, paddle"""
import deepxde as dde
import numpy as np
# 导入 deepxde 和 numpy
# 根据使用的后端不同，导入对应的深度学习框架
from deepxde.backend import tf  # 使用 tensorflow.compat.v1 作为后端
# 如果使用 paddle 后端，可以取消以下注释
# import paddle

from scipy.special import gamma  # 导入 gamma 函数，用于计算伽马函数值

# 设定分数阶微积分的阶数
alpha = 1.8

#  定义偏微分方程（FPDE）
def fpde(x, y, int_mat):
    """
    方程：\int_\theta D_\theta^\alpha u(x)
    表示在角度θ的分数阶偏微分的积分。
    
    参数：
    x - 输入位置，tensor数组，形状为（batch_size, 2）
    y - 神经网络预测的函数值，tensor（batch_size, 1）
    int_mat - 积分算子矩阵（可能是稀疏矩阵）
    """
    # 判断 int_mat 是否为稀疏矩阵，若是，则转换
    if isinstance(int_mat, (list, tuple)) and len(int_mat) == 3:
        int_mat = tf.SparseTensor(*int_mat)
        lhs = tf.sparse_tensor_dense_matmul(int_mat, y)  # 稀疏矩阵与向量乘法
    else:
        lhs = tf.matmul(int_mat, y)  # 常规矩阵乘法
    lhs = lhs[:, 0]  # 获取结果的第一列（确保shape为（batch_size,））
    
    # 计算左边项的缩放系数
    lhs *= gamma((1 - alpha) / 2) * gamma((2 + alpha) / 2) / (2 * np.pi ** 1.5)
    
    # 计算右边的已知函数值
    x = x[: tf.size(lhs)]  # 调整x的形状以匹配lhs
    rhs = (
        2 ** alpha
        * gamma(2 + alpha / 2)
        * gamma(1 + alpha / 2)
        * (1 - (1 + alpha / 2) * tf.reduce_sum(x ** 2, axis=1))
    )
    # 返回残差：左边减去右边
    return lhs - rhs

# Paddle版本的实现（已注释，若使用paddle框架可取消注释）
# def fpde(x, y, int_mat):
#     """\int_\theta D_\theta^\alpha u(x)"""
#     if isinstance(int_mat, (list, tuple)) and len(int_mat) == 3:
#         indices, values, shape = int_mat
#         int_mat = paddle.sparse.sparse_coo_tensor(
#             [[p[0] for p in indices], [p[1] for p in indices]],
#             values,
#             shape,
#             stop_gradient=False
#         )
#         lhs = paddle.sparse.matmul(int_mat, y)
#     else:
#         int_mat = paddle.to_tensor(int_mat, dde.config.real(paddle), stop_gradient=False)
#         lhs = paddle.mm(int_mat, y)
#     lhs = lhs[:, 0]
#     lhs *= gamma((1 - alpha) / 2) * gamma((2 + alpha) / 2) / (2 * np.pi ** 1.5)
#     x = x[: paddle.numel(lhs)]
#     rhs = (
#         2 ** alpha
#         * gamma(2 + alpha / 2)
#         * gamma(1 + alpha / 2)
#         * (1 - (1 + alpha / 2) * paddle.sum(x ** 2, axis=1))
#     )
#     return lhs - rhs


# 定义边界条件函数，指定在边界上 u(x) 的值
def func(x):
    # 这里的解为： (|1 - ||x||^2|) ** (1 + alpha/2)
    # ，即边界条件使用了这个表达
    return (np.abs(1 - np.linalg.norm(x, axis=1, keepdims=True))) ** (1 + alpha / 2)

# 1. 定义几何域：二维盘面，中心在[0,0]，半径为1
geom = dde.geometry.Disk([0, 0], 1)

# 2. 定义边界条件对象：在边界上 u(x) = func(x)
bc = dde.icbc.DirichletBC(geom, func, lambda _, on_boundary: on_boundary)

# 3. 构建数据对象，用于PINN训练
# 这里定义带有内部点和边界点的FPDE问题
# 参数解释：
# - 8表示内部训练点数（ull这些点在domain内部）
# - 100代表边界点数，用于强制边界条件
# - num_domain=100 表示在domain内部采样100个点
# - num_boundary=1 表示在边界采样1个点（这里用的是核心参数）
data = dde.data.FPDE(
    geom,      # 几何域
    fpde,      # 方程定义
    alpha,     # 分数阶α
    bc,        # 边界条件
    [8, 100],  # 内部点数和边界点数
    num_domain=100,  # 在domain内采样的点数
    num_boundary=1,  # 边界节点数
    solution=func  # 已知解析解（用于比较）
)

# 4. 构建神经网络
# 输入维度为2（x,y）
# 共有4层隐藏层，每层20个节点
net = dde.nn.FNN([2] + [20] * 4 + [1], "tanh", "Glorot normal")

# 5. 添加输出尺度变换（可选）
# 把网络输出乘以(1 - x^2 - y^2) ，确保边界满足特定条件（如在边界值为0）
net.apply_output_transform(
    lambda x, y: (1 - tf.reduce_sum(x ** 2, axis=1, keepdims=True)) * y
)

# 如果使用paddle，请取消下面的注释
# net.apply_output_transform(
#     lambda x, y: (1 - paddle.sum(x ** 2, axis=1, keepdim=True)) * y
# )

# 6. 将数据和网络组合成模型
model = dde.Model(data, net)

# 7. 编译模型，选择优化器和学习率
model.compile("adam", lr=1e-3)

# 8. 开始训练，训练20,000轮
losshistory, train_state = model.train(iterations=20000)

# 9. 绘制训练过程中损失变化并保存
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

# 10. 测试模型性能
# 从几何区域采样1000个点，作为测试点
X = geom.random_points(1000)
# 计算真实值
y_true = func(X)
# 用模型预测
y_pred = model.predict(X)
# 计算相对L2误差
print("L2 relative error:", dde.metrics.l2_relative_error(y_true, y_pred))

# 保存测试点的真实值和预测值到文件
np.savetxt("test.dat", np.hstack((X, y_true, y_pred)))