ARTICLE DETAIL

资讯详情

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

基于PINN的微分方程求解:PyTorch实现与避坑指南

基于PINN的微分方程求解:PyTorch实现与避坑指南 简介基于PINN的微分方程求解Python代码包面向科研人员、工程师和拥有一定Python基础的学习者系统展示物理信息神经网络求解常微分方程与偏微分问题的完整流程。内容覆盖常微分方程组、扩散方程、泊松方程、拉普拉斯方程、洛伦兹系统以及欧拉梁等典型算例同时包含DeepXDE对比实验和Jacobian-Hessian方法测试每个ipynb文件均围绕定义问题、选择网络、构建损失、训练验证展开并配有结果可视化便于从零复现并理解PINN的关键环节。压缩包共26个文件以17个ipynb示例为主体3个py脚本补充物理模型、几何域和PDE函数定义另有4个zbak备份文件、1张结果图和1份README说明整体约891KB体积小巧、结构清晰Jupyter环境即可直接运行。已有217人学习适合希望快速入门PINN并通过示例代码开展科学计算实验的中级开发者参考也可为相关课题研究提供直接复现代码。1. 基于 PINN 的微分方程求解在 Python 里落地先绕开三个误区基于 PINN 的微分方程求解方法在 Python 里落地核心就一句话让神经网络去猜解函数再用微分方程本身给这个猜测打分。第一次用的人容易踩三个误区以为它用来求解析解不是求的是数值解以为它完全免调参不是采样点、权重、激活函数都得伺候以为把方程抄进损失函数就万事大吉不是二阶导、梯度回传这些细节任一个写错训练直接虚假繁荣。这篇笔记就是照这三个误区来拆的先从损失函数讲清楚 PINN 到底在做什么再给一份能直接跑通的一阶 ODE 最小实现最后把偏微分方程场景下的参数设置、避坑记录和验证流程交代完整。适合手里有微分方程求不出解析解、又不甘心从头啃有限元剖分的工程师和科研人员。2. 把“求解”变成“优化”PINN 损失函数的最小单位拆解2.1 微分方程怎么变成残差一个一阶 ODE 的直观例子PINN 全称 Physics-Informed Neural Networks中文一般叫物理信息神经网络。这个名字已经说明白了一半神经网络负责表达解函数 u(x)微分方程本身提供监督信号。拿最简单的一阶常微分方程举例dy/dx yy(0) 1x ∈ [0, 1]解析解是 e^x没什么好算的。现在不用解析法改让一个网络来猜输入 x输出 u_θ(x)θ 是网络权重。网络是一个可微函数所以它的导数 du_θ/dx 能用自动微分精确算出来。“猜得对不对”就变成一个可以量化的指标方程残差r(x) du_θ/dx − u_θ(x)把一批 x 的采样点代进残差求平方平均就得到方程残差损失。整个训练目标就是让 r(x) 尽量趋近 0。注意这里没有拿任何“标准答案”给网络看只有方程本身的运算关系在约束它这就是“物理信息”四个字的落点方程充当了监督器。传统有限元要先剖分网格、选基函数、组装刚度矩阵解到一半想调整边界条件还得重来一遍PINN 是免网格的边界条件往损失里加一项就行几何复杂度的提升在成本上几乎无感。代价是优化过程不直观训练失败时你很难分清楚是网络容量不够、采样不合理还是权重没配平。这个“黑匣子”属性是初学者最大的挫折来源后面几章重点处理它。2.2 自动微分PINN 的导数为什么不需要手推PINN 里所有导数都不是用数值差分逼近的而是通过 PyTorch 的torch.autograd.grad按链式法则解析回传。它对 PINN 的意义在于只要网络是光滑的你要几阶导它都能给。一阶导数是这个写法du_dx torch.autograd.grad(u, x, grad_outputstorch.ones_like(u), create_graphTrue)[0]想算二阶导就在du_dx的结果上再调用一次torch.autograd.grad但第一次求导时create_graphTrue必须给否则计算图被释放第二次求导就失去了向网络权重回传的路径。这个细节是后面所有高阶导数坑的源头先记住。另外一个反直觉点你不手推导数公式不代表损失对网络参数的梯度是自动拿到的。torch.autograd.grad给出的是“输出对输入”的导数想要“损失对网络权重”的梯度还得靠常规的loss.backward()再回传一层。这条链路一旦中间哪一环断了PDE 损失就悄悄变为零训练看起来在跑实际什么都没学到。2.3 初值与边值项、数据项损失函数的完整拼图方程残差只约束区域内部初值和边值不会自动满足必须显式加项。上面一阶 ODE 的损失是这个形式L λ_pde * (1/N) * Σ r(x_i)² λ_ic * (u(0) − 1)²推广到一般问题损失函数是几块的线性组合L λ_pde * L_pde λ_ic * L_ic λ_bc * L_bc λ_data * L_data其中L_data是稀疏观测点上的拟合损失。比如某个物理过程在真实设备上只装了几个传感器把测量值作为数据项加进去PINN 就能在缺大范围数据时做参数反演这就是逆问题的入口。但项一旦多起来权重平衡就成了调试的主战场。λ_ic、λ_bc和λ_pde往往要差几个数量级而不是都填 1。具体怎么配第 4 章会给一套能直接抄的起始值。3. 用 PyTorch 跑通第一个 PINN一阶常微分方程的完整实现3.1 最小网络模型三层全连接加 Tanh 就够很多搜“PINN 代码”的人卡在第一步不知道网络该搭多深。一阶 ODE 对这种问题太简单了三层全连接、每层 50 个神经元完全够用。网络输入是 x输出是 u中间全部用 Tanh 激活。代码如下import torch import torch.nn as nn class PINN(nn.Module): def __init__(self, n_hidden3, n_neurons50): super().__init__() layers [nn.Linear(1, n_neurons)] for _ in range(n_hidden): layers.append(nn.Tanh()) layers.append(nn.Linear(n_neurons, n_neurons)) layers.append(nn.Tanh()) layers.append(nn.Linear(n_neurons, 1)) self.net nn.Sequential(*layers) def forward(self, x): return self.net(x)输入层 1 个神经元对应自变量 x输出层 1 个神经元对应 u(x)。层间宽度是n_neurons默认 50。中间夹 Tanh输出前也保留一层 Tanh这样网络输出不会因为权重初始化出现极端的尖峰值训练初期更稳。后面要解二维 PDE只需要把输入层从nn.Linear(1, n_neurons)改成nn.Linear(2, n_neurons)其余结构都不用动。这算是 PINN 最省心的地方同一个网络壳子换个损失函数就是解另一个方程。3.2 两个损失函数方程残差与初始条件的代码写法一阶 ODE 只需要两个损失函数一个是区域内的方程残差一个是初始条件。写的时候形状对齐是关键输入统一用(N, 1)的形状不要把torch.linspace出来的一维向量直接喂进去。def pde_loss(model, x): x x.clone().requires_grad_(True) u model(x) # u 对 x 求一阶导保留计算图供后续 backward du_dx torch.autograd.grad( u, x, grad_outputstorch.ones_like(u), create_graphTrue )[0] # 方程残差: dy/dx - y 0 residual du_dx - u return torch.mean(residual ** 2) def ic_loss(model, x_ic, u_ic): return torch.mean((model(x_ic) - u_ic) ** 2)这里有两处容易翻车。第一x.clone().requires_grad_(True)用于生成新的叶子张量而不是直接在采样点上原地改requires_grad避免后续训练步骤把输入错当成带梯度的网络参数。第二grad_outputstorch.ones_like(u)告诉autograd.grad对 u 的每个分量都以 1 为系数求梯度这样得到的是(N, 1)的du_dx漏掉这个参数梯度形状会退化成(N,)甚至(N, N)的雅可比初学阶段特别容易在这里报形状错误。注意create_graphTrue在这两个损失里必须写。PDE 损失要经过残差回到网络权重如果这里把计算图丢了loss.backward()要么报错要么 Pde 项梯度为空看起来训练在跑实际上只有初始条件在更新。3.3 训练循环与权重设置为什么初始条件要乘 10采样点和损失项都准备好之后训练循环本身很朴素。下面这份代码把方程残差和初值条件组合在一起lambda_ic 10.0这个权重不是拍脑袋背后是点数量级的差异。model PINN() optimizer torch.optim.Adam(model.parameters(), lr1e-3) # 教学场景用固定采样点实际项目建议每轮重新采样 x_pde torch.linspace(0, 1, 300).view(-1, 1) x_ic torch.zeros(1, 1) u_ic torch.ones(1, 1) lambda_ic 10.0 for step in range(15000): optimizer.zero_grad() loss pde_loss(model, x_pde) lambda_ic * ic_loss(model, x_ic, u_ic) loss.backward() optimizer.step() if step % 1000 0: print(fstep {step:5d} loss {loss.item():.3e}) # 检查初值有没有被压住 print(u(0) , model(x_ic).item())方程残差有 300 个采样点初始条件只有 1 个点。如果不加权pde_loss的梯度比ic_loss大两个数量级Adam 会优先压残差初值误差很容易停在 0.1 这个级别。把lambda_ic提到 10 甚至 50初始条件才压得住。后面换成 PDE 场景边界条件也同样需要这个思路不然你会看到 loss 降得很漂亮但在端点一验证边界值就是不对。诊断技巧每 1000 步打印 loss 的同时打印u(0)。如果 loss 快速下降但u(0)离 1 越来越远就是lambda_ic太小直接往上加就行。训练到 15000 步时这个模型在 [0,1] 内的最大误差通常能到 1e-5 以下对一个入门算例来说已经够用。4. 从 ODE 扩到偏微分方程网络、激活、采样与优化器的四个关键改动4.1 以 Poisson 方程为例二维输入网络与二阶导写法ODE 跑通之后扩到 PDE 不需要换框架只需要改四处输入维度、导数阶数、边界项、采样密度。以二维 Poisson 方程为例−u_xx − u_yy f(x, y)定义在 [0,1]²边界 u 0网络输入从一维变成两个坐标的拼接残差则要对 x 和 y 各求二阶导。损失函数写法如下def pde_loss_2d(model, xy, f_source): x xy[:, 0:1].clone().requires_grad_(True) y xy[:, 1:2].clone().requires_grad_(True) u model(torch.cat([x, y], dim1)) u_x torch.autograd.grad(u, x, grad_outputstorch.ones_like(u), create_graphTrue)[0] u_y torch.autograd.grad(u, y, 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] u_yy torch.autograd.grad(u_y, y, grad_outputstorch.ones_like(u_y), create_graphTrue)[0] residual -u_xx - u_yy - f_source(x, y) return torch.mean(residual ** 2)两次调用autograd.grad都必须写create_graphTrue。第一次求u_x时开图是为了第二次能继续求u_xx第二次求u_xx时开图是为了残差能反向传播到网络权重。任何一次漏掉二阶导相关的那一项梯度都会静默丢失。f_source是方程右端项测试时常用常数 f1能看出解在中心处下凹。边界损失单独写一个函数把四条边界的点集拉出来计算model(xy_bc)²的均值。边界点数量建议和内部点在同一量级别出现内部点 10000 个、边界点只有 100 个的情况那样边界条件会被稀释。4.2 激活函数为什么不能随便换ReLU 在二阶方程里的陷阱激活函数对 PINN 的影响远大于普通监督学习。Tanh 是默认选择因为它光滑、二阶导连续、输出范围有限正好匹配微分方程对光滑性的要求。sin 在带周期性的问题上效果不错但对权重初始化和学习率敏感容易振荡。ReLU 在这里要重点提醒它分片线性二阶导恒等于 0。放进 Poisson 损失里u_xx和u_yy永远计算为 0残差只剩下-f和网络输出毫无关系方程信息直接丢失。你可能会看到边界损失在下降但内部完全是一团乱。真要用 ReLU 系的激活只能用平滑版本比如 SiLU 或 GELU。我的经验是二阶 PDE 场景别折腾直接 Tanh最多把网络加深到 4 到 6 层去换表达力而不是换一个听起来更高级的激活函数。4.3 采样策略固定网格、随机重采样与残差自适应加点采样密度直接决定 PINN 的解质量。一阶 ODE 用linspace固定 300 个点没问题二维 PDE 强烈建议每次迭代重新随机采样。做法是每个 step 从 [0,1]² 均匀生成一批新点网络见到的是整个分布而不是某组固定坐标。固定坐标有个很隐蔽的坑网络本质上在拟合一组确定坐标训练点之间的行为约束很弱验证时只要取样点偏一点误差立刻暴露。你看到训练 loss 到 1e-6以为收敛了换一网格点一算PDE 残差可能还在 1e-2 量级。更强力的是残差自适应加点每 500 步在当前解上算每个点的残差绝对值找出 top 10% 的大残差点以它们为中心在半径 ε 内补一批随机点下一轮混进训练集。这个策略把网络容量集中在方程最难满足的区域激波、边界层这类问题会明显受益。4.4 优化器顺序Adam 热身加 LBFGS 精修的真实用法PINN 优化器的主流玩法是两段式Adam 先把损失压到平台期再切 LBFGS 精修。Adam 在前期不容易炸适合大范围搜索LBFGS 是拟牛顿法对光滑残差面的收敛更准经常能把损失从 1e-4 推到 1e-6 甚至更低。但 LBFGS 有两条硬约束必须是全批量计算不能配合 mini-batch 随机采样学习率要保守从 0.1 起步。切换代码optimizer_lbfgs torch.optim.LBFGS(model.parameters(), lr0.1, max_iter20) def closure(): optimizer_lbfgs.zero_grad() loss pde_loss(model, x_pde) lambda_ic * ic_loss(model, x_ic, u_ic) loss.backward() return loss for _ in range(50): optimizer_lbfgs.step(closure)如果 LBFGS 一进去 loss 就反弹把 lr 从 0.1 降到 0.01多数情况下能稳住。注意 LBFGS 内部步数未必等于我们的外层循环次数观察指标是外层每次epoch后 loss 是否单调下降而不是看迭代次数。参数设置可以直接参照这张表问题类型网络结构采样策略优化器一阶 ODE3 层 × 50Tanh200~500 点固定或重采样Adam lr1e-3二阶 ODE / 低维稳态 PDE4 层 × 80Tanh1000~3000 点重采样Adam 热身 LBFGS含时间 PDE / 强非线性5~6 层 × 100Tanh 或 sin5000 点残差自适应Adam 热身 LBFGS含时间的 PDE 不需要特殊的网络结构把时间 t 当普通输入特征和空间坐标拼在一起卷积都不用加损失函数按原方程写就行。5. PINN 训练避坑与排查五个让我翻车的真实原因PINN 训练不像普通监督学习那么直观网络翻车时你很难判断是哪一环出了问题。下面五条是按出现频率排的每条都按现象、原因、解决的顺序写都是我实际踩过且找到明确根源的问题。5.1 PDE 损失不下降只有初值项在动现象训练几千步pde_loss纹丝不动只有ic_loss在下降或者 loss 整体下降但解函数完全不符合方程形态。原因最常见的是激活函数用了 ReLU二阶导恒为零PDE 残差变成一个常数梯度算不出来。另一个高频原因是torch.autograd.grad里漏了create_graphTrue一阶导以上的计算图断开PDE 项反向传播不了。解决换成 Tanh 激活并且所有autograd.grad调用都显式写create_graphTrue。排查技巧是打印loss_pde.item()和loss_ic.item()两个分量而不是只看总 loss一眼就能看出哪项在偷懒。5.2 训练 loss 很小换一组点验证却对不上现象训练时采样点上都对loss 压到 1e-6把验证点取在训练点之间或区间边缘误差突然涨到 1e-2。原因训练点固定且分布规则时网络本质上在记忆这组坐标。两个训练点之间的行为只有方程残差在约束残差对每个点贡献不均匀网络会在采样稀疏区域“偷懒”。解决每个 step 重新随机采样别用固定网格。验证时也要用全新点算残差不要看训练 lossx_new torch.rand(1000, 1) # 全新随机点 loss_on_new pde_loss(model, x_new).item() print(loss_on_new)如果loss_on_new比训练 loss 大一个量级以上说明采样不够或网络过拟合了固定点优先补采样。5.3 边界条件在训练尾声悄悄变差现象前几千步初值和边值都很准训练上万步之后u(0)从 1.0 漂到了 0.98 左右且这个漂移不反弹。原因PDE 残差采样点数量远超初值/边值点训练后期学习率变小梯度被残差项主导边界点被稀释。每轮更新的重心都在内部区域边界成了盲区。解决把lambda_ic提到 50~100或者把初值点从单个点改成初值附近的小邻域比如在 [-0.01, 0.01] 内取 5 个点梯度更稳。更稳妥的做法是在每 500 步训练后单独打印一次model(x_ic)盯住边界误差别等训练结束再惊讶。5.4 输出在采样点之间剧烈振荡现象loss 正常下降但把预测的 u(x) 画出来曲线呈波浪形或锯齿形导数一高一低明显不光滑。原因网络容量偏大、训练点过少或者激活函数用了 sin 且学习率偏高。LBFGS 的学习率太大也会出现这种情况一步迈过头在陡峭区域来回震荡。解决先把网络宽度从 100 减回 50增加采样点密度激活换回 TanhLBFGS 的 lr 收到 0.01 试一遍。原则上先减容量而不是加正则PINN 的振荡大多是容量和采样密度不匹配导致的。5.5 二阶导模型的梯度爆炸与 NaN现象算 Poisson 这类二阶方程时loss 越训越大甚至直接变成 NaN。原因二阶导数对网络权重的敏感度远高于一阶损失曲面更陡。学习率给到 1e-2 基本必炸网络初始化方差偏大时Tanh 的饱和区和二阶导叠加梯度会失控。解决学习率降到 1e-4 到 1e-3网络权重用nn.init.xavier_normal_初始化把 x 从 [0,1] 归一化到 [-1,1]能显著缓解高阶导数的数值问题。归一化这步很容易被忽略但它对二阶 PDE 的帮助几乎是立竿见影的。6. 验证一套 PINN 结果可不可信从解析解对齐到数值差分交叉验证6.1 先用带解析解的方程把全流程校准一遍PINN 里“loss 低”不等于“解对了”。我拿到一个新方程时的第一件事不是直接在目标方程上调参而是找一个同类型、带解析解的方程把整个流程跑通。比如后面要解的是二维稳态热传导就先用一个已知温度场的 Poisson 方程做基准。验证代码很简单把训练域外扩一点测试x_test torch.linspace(-0.5, 1.5, 1001).view(-1, 1) u_pred model(x_test).detach().numpy().ravel() u_true np.exp(x_test.numpy().ravel()) # dy/dx y 的解析解 print(max abs err:, np.max(np.abs(u_pred - u_true)))测试范围比训练区间宽PINN 的外插能力有限这个误差通常会比区间内大不少。但这一步能告诉你训练的边界在哪如果连训练区间内部都对不齐到 1e-6问题一定出在训练环节而不是验证环节。6.2 在新采样点上算残差而不是看训练 loss验证残差要用新随机点。做法和训练时一样只是这次不更新梯度纯粹检查方程被满足的程度x_check torch.rand(2000, 1) * 1.2 - 0.1 u model(x_check) du torch.autograd.grad(u, x_check, grad_outputstorch.ones_like(u), create_graphTrue)[0] res (du - u).detach().numpy() print(rmse:, np.sqrt(np.mean(res**2)), p99:, np.percentile(np.abs(res), 99))只看 rmse 不够要看 p99 分位数。残差分布的尾部才是问题点集中区。如果 p99 明显大于 rmse说明大部分点都满足方程但有一小块区域严重不满足优先去那里补采样点。网格类方法看残差最大值PINN 必须看残差分布这是两者思维上最大的差别。6.3 验证完再做优化我现在的固定习惯我现在的习惯是PINN 训练完之后先做两件笨事一看初值和边值偏差二在新采样点上做残差直方图。这两件事做完才敢把结果交给下游计算。之前有一次把训练 loss 压到 1e-7 的“漂亮模型”放进对比基准结果一查边值偏了 0.2整个基准作废。从那以后验证永远排在调参前面。数值方法本身的病态问题往往藏在你看不见的间隙里先证明这套流程在简单问题上可靠再拿去解复杂问题才是让 PINN 从“玄学”变成工具的唯一路径。希望帮到你。本文还有配套的精品资源点击获取
返回列表