ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

PINN物理信息神经网络入门:用PyTorch求解微分方程的完整实战

PINN物理信息神经网络入门:用PyTorch求解微分方程的完整实战 简介面向科学计算与深度学习交叉领域的 Python PINN 求解微分方程资料包适合具备一定 Python 与深度学习基础希望用物理信息神经网络替代传统数值方法在科研或工程中快速验证微分方程解的工程师、研究生与科研人员。资源共 26 个文件以 17 个可运行的 ipynb 案例笔记本为主py 脚本用于定义网络结构与损失函数README 梳理方法流程效果图可辅助对比不同边界条件下的求解表现压缩包约 891KB。目录覆盖 Euler 梁、扩散方程、Laplace 方程、Poisson 方程、常微分方程系统、Lorenz 系统等典型问题并包含 Jacobian-Hessian 方法测试与 DeepXDE 库实践展示不同实现路线。每个案例笔记本均包含从数据准备、网络结构设计、损失函数构建到训练与结果分析的完整流程且针对不同方程给出 Dirichlet、Neumann、Robin、Periodic 等边界条件处理、导数约束方式与优化器选择可直接运行并改写迁移到自己的方程场景。已有 217 人学习适合作为 PINN 入门到进阶的实践参照帮助快速掌握自动微分与物理约束建模的关键思路。1. 用神经网络解微分方程PINN为什么值得学提到微分方程大多数人第一反应是有限差分、有限元这些老牌数值方法它们要画网格、组装矩阵、迭代求解换个边界条件可能整套流程重来。PINNPhysics-Informed Neural Network物理信息神经网络的思路完全不同它把微分方程本身当作训练误差的一部分让神经网络在拟合初边值条件的同时尽量满足方程残差为零。你不需要剖分网格只需要在求解域内撒点就能得到一个连续、可求导的近似解。它尤其适合复杂几何、逆问题、稀数据场景。对做Python工程的人来说PINN的门槛不高会写PyTorch基本训练循环就能很快跑通第一个模型。这篇笔记我会从损失函数构造讲到最小可跑代码再讲训练中的玄学和踩坑最后落到验证手法上。2. PINN的数学框架把微分方程变成一种“可训练的错误”2.1 损失函数不是一项残差、边界、初始、数据各司其职传统神经网络训练分类任务损失函数只有一项比如交叉熵。PINN却要同时伺候多个“老板”。以一个一般的含时间项微分方程为例假设方程为[ \mathcal{L}[u(x)] f(x) ]其中 \mathcal{L} 是对 u 的某种微分算子。你用一个神经网络 \hat{u}(x; \theta) 去逼近真实解 u那么我们希望这个 \hat{u} 满足两个东西第一带进去方程左边等于右边即残差 \mathcal{L}[\hat{u}] - f(x) 尽可能接近 0第二它满足问题的初始条件和边界条件。所以训练时分母上的损失至少包含三块残差损失PDE Loss在内部采样点 x_i 上计算 \frac{1}{N_f}\sum_{i}\left(\mathcal{L}[\hat{u}(x_i)] - f(x_i)\right)^2。边界损失BC Loss在边界点 x_b 上计算 \left(\hat{u}(x_b) - g(x_b)\right)^2。初始损失IC Loss瞬态问题里在 t0 时刻采样计算 \left(\hat{u}(x,0) - h(x)\right)^2。如果还有观测数据比如某个位置的真实测量值 u_{obs}再加上数据损失项 \left(\hat{u}(x_{obs}) - u_{obs}\right)^2。这里最容易犯的错是把三项做成一个大 loss 然后平均因为各区域采样点数量不平衡时项与项会有天然的尺度差。边界点只有几十个内部点可能上千个如果直接求和内部残差主导训练边界条件可能被彻底淹没。我一般先把每一项归一化再乘以权重后面第四章细说。2.2 自动微分PINN 的“求导引擎”和普通神经网络的区别普通神经网络训练反向传播要计算损失对网络参数的梯度。PINN 额外需要的是网络输出 \hat{u} 对输入坐标 x 的偏导数而且可能是一阶、二阶甚至高阶。这些导数通过链式法则一层一层传回去正是 PyTorch 和 TensorFlow 自动微分autograd引擎擅长的。举个例子计算 \hat{u} 对 x 的一阶偏导在 PyTorch 里就是u model(x) u_x torch.autograd.grad(u, x, grad_outputstorch.ones_like(u), create_graphTrue, retain_graphTrue)[0]参数解释u网络输出shape 为(N,1)。x输入坐标shape 为(N,1)要求 requires_gradTrue否则无法求导。grad_outputs因为u是向量torch.autograd.grad需要给一个和u同形状的初始梯度向量一般填全 1表示对每个输出分量各自求梯度。create_graphTrue必须打开因为这个导数后面还要参与二次求导如二阶导数并且要用于损失对参数的梯度计算需要把计算图保留下来。retain_graphTrue多次调用grad时避免计算图被释放否则第二次求导会报错。更高阶导数就是在这个一阶导数上再做一次同样的操作。这个“求导引擎”直接决定了 PINN 能不能跑起来。如果你用数值差分去算这些导数不仅慢而且误差会污染损失导致训练永远没法收敛到机器精度。2.3 损失权重四个目标打架时怎么拉架PINN 训练中常遇到一个现象残差损失已经降到 1e-5但边界损失还在 1e-2总体 loss 看起来还是很大。或者反过来边界拟合得很好但方程内部完全不对。这就是“多目标打架”。残差、初边值、数据这几项它们的梯度尺度可能差好几个数量级。比如一个变化剧烈的解它的二阶导数值很大残差 loss 的梯度也大就会把网络参数往“满足方程但破坏边界”的方向猛推。常见做法是给每项加一个权重系数\text{loss} \lambda_{pde} L_{pde} \lambda_{bc} L_{bc} \lambda_{ic} L_{ic}。但这些 \lambda 怎么设全靠经验和试。我个人的起点值是\lambda_{pde}1.0\lambda_{bc}10 \sim 100\lambda_{ic}10 \sim 100。因为边界点远少于内部点必须放大边界权重。后续如果发现总体 loss 明明在降但边界误差不减就把边界权重再向上调 10 倍。更现代一点的做法是“自适应权重”每训练若干个 step重新统计各项损失的梯度范数然后按梯度范数成反比分配权重这就是 GradNorm 思路。它可以省去手动调参的苦力但也带来额外的计算开销。对于刚接触 PINN 的人我建议先从固定权重开始跑通一个简单问题之后再上自适应。3. 用 PyTorch 跑通第一个 PINN以一阶 ODE 为例3.1 环境准备与依赖torch、numpy、matplotlib做 PINN 只需要基础科学计算三件套。PyTorch 负责网络定义和自动微分NumPy 做数据生成Matplotlib 画验证曲线。安装时注意 Python 版本和 CUDA 版本匹配如果只是跑 CPU 版本直接pip install torch --index-url https://download.pytorch.org/whl/cpu可以避开庞大的 CUDA 依赖。我建议新手先在 CPU 上把流程走通因为一维 ODE 的算力需求很小CPU 几秒就能完成一次训练没必要一上来就折腾 GPU 环境。另外要注意 PyTorch 和 Python 的版本兼容性。比如 Python 3.12 刚出来时某些老版本 PyTorch 没有对应 wheel安装会直接报错。如果你用 Anaconda 管理环境可以用conda create -n pinn python3.10这个版本兼容性最稳。3.2 完整代码定义网络、损失、训练循环我们求解一个最简单的常微分方程[ \frac{dy}{dt} -y, \quad y(0)1 ]解析解是 y e^{-t}方便验证。PINN 要做的事在 t 属于 [0, 4] 区间内采样一堆点让网络输出 y(t) 满足方程导数关系并且满足初始条件。下面是一份完整可跑的最小代码。import torch import torch.nn as nn import numpy as np import matplotlib.pyplot as plt # 1. 定义网络输入一个数(t)输出一个数(y) class PINN(nn.Module): def __init__(self): super(PINN, self).__init__() self.net nn.Sequential( nn.Linear(1, 20), nn.Tanh(), nn.Linear(20, 20), nn.Tanh(), nn.Linear(20, 20), nn.Tanh(), nn.Linear(20, 1) ) def forward(self, t): return self.net(t) # 2. 定义损失 def pde_loss(model, t): t.requires_grad_(True) y model(t) y_t torch.autograd.grad(y, t, grad_outputstorch.ones_like(y), create_graphTrue)[0] # 方程 dy/dt -y 残差 y_t y residual y_t y return torch.mean(residual**2) def ic_loss(model, t0): y0 model(t0) return torch.mean((y0 - 1.0)**2) # 3. 训练准备 model PINN() optimizer torch.optim.Adam(model.parameters(), lr1e-3) lambda_pde 1.0 lambda_ic 10.0 # 内部采样点区间 [0, 4] 均匀取 200 个点 t_f torch.linspace(0, 4, 200).reshape(-1, 1) # 初始条件点只用 t0 t_0 torch.tensor([[0.0]]) # 4. 训练循环 loss_history [] for epoch in range(5000): optimizer.zero_grad() loss_pde pde_loss(model, t_f) loss_ic ic_loss(model, t_0) loss lambda_pde * loss_pde lambda_ic * loss_ic loss.backward() optimizer.step() loss_history.append(loss.item()) if epoch % 500 0: print(fepoch {epoch}, pde_loss: {loss_pde.item():.6f}, ic_loss: {loss_ic.item():.6f}) # 5. 验证 t_test torch.linspace(0, 4, 100).reshape(-1, 1) y_pred model(t_test).detach().numpy() y_true np.exp(-t_test.numpy()) plt.plot(t_test.numpy(), y_true, labelexact) plt.plot(t_test.numpy(), y_pred, --, labelPINN) plt.legend() plt.show()逻辑说明在pde_loss里我让输入t开启requires_grad然后通过网络输出y再求y对t的一阶导数y_t。这样residual y_t y就对应原方程移项后的左边。注意这里方程是 dy/dt -y所以残差写成 y_t y 等于 0。如果你换一个方程残差表达式要对应改。在训练循环里每一项损失独立计算加权求和后再反向传播。内部采样点用的是固定t_f初始条件点就是一个单独的t_0。固定采样点简单但会在训练中过拟合这些点后面我会讲“每轮重新采样”的变体。初始条件权重lambda_ic10是我试出来的。如果设为 1前几百步网络会先费劲地满足内部方程初始条件反而滞后调大后训练初期网络会优先把 y(0) 压到 1再逐步修正内部。实际训练几千步后 pde_loss 会降到 1e-4 以下ic_loss 也会降到 1e-5 左右。3.3 训练完怎么判断“解对了”三个验证指标损失降下来不代表解就对了。PINN 有一个特点损失函数只是方程和边界的平均平方误差它可能因为采样点稀疏而在未采样位置悄悄“作弊”。所以验证必须用独立的稠密样本。第一个指标是解析解对比。比如上面这个 ODE把 t 加密到 1000 个点计算最大绝对误差 \max |y_{pred} - e^{-t}|。如果这个值小于 1e-3基本算成功。第二个指标是残差在验证点上的分布。重新随机采 500 个点跑一遍 pde_loss看看是不是处在同一数量级。如果验证集残差远大于训练集残差说明过拟合到了固定采样点。第三个指标是导数一致性。对瞬态问题把网络输出的导数曲线和解析导数曲线画在一起观察是否贴合。有时候 y 本身拟合得好导数却差得远说明网络在“硬凑”函数值没有真正学到方程结构。4. 训练策略与参数调优PINN 的“手感”从哪来4.1 网络结构深度宽度和激活函数的取舍PINN 的标配激活函数是 Tanh不是 ReLU。原因很简单微分方程需要有足够光滑的导数ReLU 在 x0 处不可导其导数分段常数二次导数为零二阶导数的信息全丢无法表达曲率丰富的解。而 Tanh 无穷阶光滑适合做数值微分。你如果用 ReLU 去解扩散方程含 u_{xx}残差会大面积变成 0训练毫无进展。网络结构上我一般从 3 层隐藏层、每层 20 个神经元起步。这个配置能解大多数一维基准问题。如果解很光滑可以减到 2 层 16 个神经元如果解有陡峭峰或边界层要加到 4~6 层、每层 40~60 个神经元。加深比加宽更有效因为每层非线性组合能构造更复杂的导数结构但注意层数过深会带来梯度消失训练更不稳定。层数和宽度的选择在你的具体问题上没有万能公式。一个实用原则先跑一个中等规模网络看验证误差如果验证误差降不下去加宽如果训练不稳定、loss 剧烈震荡加深但加小的残差连接ResNet 风格会更稳。PINN 问题通常没有海量数据集所以不用担心过拟合数据集主要担心表达力不够。4.2 采样点与学习率先粗后细还是固定采样采样点决定了网络在哪些位置“看”到方程。固定采样上面代码那种简单但容易造成局部欠采样。更稳的方式是每轮迭代重新随机采样称为“逐轮重采样”。这样做可以避免网络记住固定点还能让每个区域被均匀覆盖的概率更高。实践中我会用一个折中策略训练前 1000 轮用固定采样让网络快速抓住全局形状之后每 500 轮重新采样一次让网络去修正细部。采样点数可以这样估一维问题 200~500 点二维问题每轴 20~30 个网格展开再加 1000 个随机点三维问题用 Halton 序列或拉丁超立方采样保证均匀性。不要用简单均匀网格尤其在高维问题上均匀网格的“对角线”区域采样稀疏会导致解在角落失真。学习率方面Adam 的默认 1e-3 在 PINN 上是比较安全的起点。如果训练 loss 下降太慢可以把学习率调到 2e-3如果出现振荡或 NaN调到 5e-4。更进阶的做法是学习率衰减每 1000 轮乘以 0.5让网络后期做细粒度收敛。这种分阶段策略比一直用固定学习率保底效果好。4.3 损失权重动态调整从固定权重到梯度归一化前面讲到固定权重 \lambda 是个“拉架”的角色但 \lambda 一旦固定就会出现“前期合适、后期不合适”的问题。训练初期残差损失往往很大比如 10边界损失较小比如 0.1此时边界权重就算给 10总占比也不高训练后期反过来残差损失变得很小边界权重 10 会让边界项占比过大训练被边界主导内部方程又开始反弹。解决这个问题最简单的办法是让权重跟着损失尺度走即每一轮计算完各项损失后把权重调整为某项损失倒数的开方或者做梯度归一化。这里给出一个我常用的“梯度归一化”核心片段# 假设 loss_pde, loss_ic 已经计算出来 # 计算各自对模型参数的梯度范数 grad_pde torch.autograd.grad(loss_pde, model.parameters(), retain_graphTrue, create_graphFalse) grad_ic torch.autograd.grad(loss_ic, model.parameters(), retain_graphTrue, create_graphFalse) norm_pde torch.sqrt(sum([g.norm()**2 for g in [_g for _g in grad_pde if _g is not None]])) norm_ic torch.sqrt(sum([g.norm()**2 for g in [_g for _g in grad_ic if _g is not None]])) # 权重与梯度范数成反比再加一个缩放系数 lambda_pde 1.0 / (norm_pde 1e-6) lambda_ic 1.0 / (norm_ic 1e-6) # 也可以进一步归一化到和为1 total lambda_pde lambda_ic lambda_pde, lambda_ic lambda_pde / total, lambda_ic / total注意torch.autograd.grad计算梯度时会额外产生一次反向传播开销但换来的是权重的自动平衡。这里retain_graphTrue是必须的因为后面loss.backward()还要继续用计算图。求梯度范数时把每层参数的梯度范数平方和开根号忽略 None 梯度比如未参与该损失的参数。使用这种策略后训练的前期和后期权重都能自行调节整体 loss 的曲线更平稳。代价是每轮训练耗时多了 20~30%在物理上可接受。如果你的问题本身很简单固定权重足够没必要上这套。4.4 训练终止条件什么时候可以停PINN 的 loss 和传统机器学习 loss 不太一样降到多低算“解对了”没有统一标准。我习惯设定一个组合条件pde_loss、bc_loss、ic_loss 各自低于某个阈值比如 1e-5并且验证集的最大绝对误差不再明显下降。记录验证误差每 200 轮的变化如果连续 5 次验证误差变化小于 1%就提前停止。不要只看总 loss 的绝对值因为不同权重下总 loss 没有可比性。训练轮数的经验值一维 ODE 3000~5000 轮足够二维稳态 PDE 一万轮左右瞬态问题几万轮起步。如果加了自适应采样轮数需求会上升因为每轮重新采样会让网络不断面对新的数据点收敛会慢但最终精度更高。用早停法时记得保存最佳模型权重别让最后一步把最优解震荡掉了。5. PINN 常见问题与避坑我踩过的四个坑5.1 训练不收敛损失停在某个平台期现象loss 在前几百轮下降后卡在 1e-2 左右怎么训练都纹丝不动。原因最常见是初始条件和边界条件没被“看见”。如果边界权重太小网络优先满足内部残差边界点上的误差始终保留但因为边界点数量少对梯度贡献小所以整体 loss 停在平台。还有一种可能是采样点分布太差比如边界点只采样了一个而内部点密度极高。解决先检查各项损失的数值。把每项单独打印出来看是哪个项卡住。如果是边界项卡住把对应权重调大 100 倍试试。如果是残差项卡住检查残差表达式是否写错比如符号错了、导数求少了。在代码里写个简单测试用解析解代入 pde_loss如果输出不是接近 0说明损失函数定义有 bug。5.2 解在边界漂亮内部却发疯现象边界条件拟合得很好但解在内部区域出现尖锐的毛刺、振荡或者完全偏离真实解。原因内部采样点数不足或者激活函数选择不当。Tanh 虽然光滑但网络深度不够时很难同时满足边界和内部复杂曲率。另外训练初期网络参数随机此时边界附近的局部梯度非常大会把网络“吸”向边界拟合内部区域长时间得不到关注。解决把内部采样点从 200 增加到 1000并使用逐轮重采样。同时把学习率调低到 5e-4让网络更稳定地探索。还有一个窍门先把内部采样点分成两个集合一固定一随机固定集用于稳定全局形状随机集用于平滑局部细节。我屡试不爽。5.3 梯度爆炸或出现 NaN现象训练到一半 loss 突然变成 nan或者梯度绝对值超过 1e10。原因学习率过大、网络初始化不当、或损失函数里存在除零操作。PINN 的损失包含高阶导数如果方程本身存在奇异性比如解在某个点导数无穷大梯度也会爆炸。另外Adam 自适应学习率虽然能缓解但遇上极端梯度依旧溃败。解决采用 Xavier/He 初始化PyTorch 默认但更稳的是把网络输入归一化到 [-1,1] 或 [0,1]。如果你求解域是 [0, 100]直接送进去初始梯度可能大得离谱。把 t 缩放成 t/100损失函数里也要相应调整导数链式法则dy/d(t_scaled) 100 * dy/dt_scaled。其次把学习率降到 1e-4梯度裁剪torch.nn.utils.clip_grad_norm_(model.parameters(), max_norm10.0)这行代码放在loss.backward()之后、optimizer.step()之前。它会把所有参数的梯度的总范数剪裁到不超过 10能有效防止个别的极端梯度带偏方向。5.4 结果随随机种子漂移换台机器数字就变现象同一个代码两次运行的验证误差相差 10 倍甚至解的形状都不一样。原因网络初始化随机、采样点随机PINN 的损失面存在大量局部最优。初值和随机点不同会收敛到不同的“坑”里。解决固定随机种子是一个办法但更实用的是“多种子集成”。跑 5 次不同的随机种子对最后的结果取平均或者取误差最小的那一次。平均可以消除部分随机振荡取最好则能给你“这个架构上限在哪”的判断。如果你追求可复现性在代码开头设置torch.manual_seed(42) np.random.seed(42)但这只保证同一台机器、同一环境可复现。跨机器仍然可能因浮点运算顺序不同产生细微差异。所以我在工程上从不用单次结果下结论至少跑三次。5.5 数据点坐标范围差太多PINN 直接罢工现象二维问题中 x 范围是 [0,1]y 范围是 [0,1000]PINN 训练不收敛或收敛后误差主要集中在大尺度维度。原因神经网络对输入特征的尺度很敏感。两个维度放一起y 尺度是 x 的千倍网络初始输出被大尺度坐标主导小尺度维度的信息几乎淹没。物理方程的 term 之间也会因为尺度不平衡导致残差计算失真。解决在进入网络之前把每个维度归一化到 [-1,1] 或 [0,1]。输出层看情况如果解的物理量级很大可以在输出后乘一个常数。注意如果你做了输入缩放自动微分求导时要把缩放因子补偿回来。例如 t_scaled t / T那么 \frac{d u}{d t} \frac{d u}{d t_scaled} \cdot \frac{1}{T}。这个补偿若是漏了你会发现网络拟合的函数在物理尺度上完全不对损失却很小。6. 从复现到实用让 PINN 在更真实的问题上站稳6.1 给网络加一个时间维度瞬态问题自然解瞬态问题比如热传导方程 \frac{\partial u}{\partial t} - \alpha \frac{\partial^2 u}{\partial x^2} 0网络输入变成 (x, t) 两个变量输出仍然是 u。损失函数里残差项要同时计算时间一阶导和空间二阶导def pde_loss(model, x, t): u model(x, t) u_t torch.autograd.grad(u, t, grad_outputstorch.ones_like(u), create_graphTrue)[0] u_x torch.autograd.grad(u, x, grad_outputstorch.ones_like(u), create_graphTrue)[0] u_xx torch.autograd.grad(u_x, x, grad_outputstorch.ones_like(u_x), create_graphTrue)[0] residual u_t - alpha * u_xx return torch.mean(residual**2)这里有一个隐藏细节u_x对应的是网络输出对输入 x 的偏导它本身也是一个张量可以作为后续grad的输入但注意create_graphTrue必须保留。采样时x 在空间域采样t 在时间域采样同时还要加上初始条件t0和边界条件x边界的采样。6.2 自适应采样服务陡峭解和局部奇异性均匀采样在高梯度区域分辨率不足尤其当解有锐利锋面或边界层比如流体力学中的激波、对流扩散方程中的小扩散系数情况。这时候要引入残差自适应采样训练一段时间后在现有模型的残差较大区域加密采样点。做法每 1000 轮跑一次当前模型在候选点上的残差绝对值。按残差大小或残差绝对值乘以某个幂作为概率分布重新采样下一批训练点。常见候选点从均匀分布中抽取数量是训练点数的 5~10 倍。代码示意with torch.no_grad(): # t_candidate 有 5000 个均匀候选点 residual_abs torch.abs(y_t y) # 每个候选点的残差绝对值 prob residual_abs / residual_abs.sum() idx torch.multinomial(prob, num_samples200, replacementTrue) t_new t_candidate[idx]注意torch.multinomial输入的prob需要一维浮点张量且所有元素非负、和为 1。这样采样频率直接和残差大小挂钩相当于让网络把“力气”花在最需要修正的地方。这种自适应策略在解存在奇异性的问题上非常有用能把边界层误差压低两三个数量级。6.3 用误差图和收敛曲线做最终验收最后验收不只看 loss我每次会画两张图。第一张是验证集上的绝对误差空间分布比如 \log_{10}|u_{pred} - u_{ref}| 的等高线图哪里亮就说明哪里误差大。第二张是训练过程中验证误差和各项损失的曲线横轴是 epoch纵轴是对数尺度。如果曲线在最后阶段还明显下降说明训练不足如果曲线平稳但验证误差不如预期说明模型表达力或采样策略有问题而不是再训就能解决。做完这一切我养成的习惯是先拿一个带解析解的简单问题验证整个流程再上真实工程问题。因为 PINN 的隐藏坑很多但如果你在简单问题上能稳定达到 1e-4 量级误差那么复杂问题上的失败至少可以排除代码 bug 和参数基础设置的问题。希望我这套“先搭框架、再调手艺、最后拿误差剖面验收”的流程能帮你在自己的微分方程求解任务里少走几趟弯路。本文还有配套的精品资源点击获取
返回列表