简介:基于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 = y,y(0) = 1,x ∈ [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_outputs=torch.ones_like(u), create_graph=True)[0]想算二阶导,就在du_dx的结果上再调用一次torch.autograd.grad,但第一次求导时create_graph=True必须给,否则计算图被释放,第二次求导就失去了向网络权重回传的路径。这个细节是后面所有高阶导数坑的源头,先记住。
另外一个反直觉点:你不手推导数公式,不代表损失对网络参数的梯度是自动拿到的。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_hidden=3, n_neurons=50): 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_outputs=torch.ones_like(u), create_graph=True )[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_outputs=torch.ones_like(u)告诉autograd.grad对 u 的每个分量都以 1 为系数求梯度,这样得到的是(N, 1)的du_dx;漏掉这个参数,梯度形状会退化成(N,)甚至(N, N)的雅可比,初学阶段特别容易在这里报形状错误。
注意:
create_graph=True在这两个损失里必须写。PDE 损失要经过残差回到网络权重,如果这里把计算图丢了,loss.backward()要么报错,要么 Pde 项梯度为空,看起来训练在跑,实际上只有初始条件在更新。
3.3 训练循环与权重设置:为什么初始条件要乘 10
采样点和损失项都准备好之后,训练循环本身很朴素。下面这份代码把方程残差和初值条件组合在一起,lambda_ic = 10.0这个权重不是拍脑袋,背后是点数量级的差异。
model = PINN() optimizer = torch.optim.Adam(model.parameters(), lr=1e-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(f"step {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], dim=1)) u_x = torch.autograd.grad(u, x, grad_outputs=torch.ones_like(u), create_graph=True)[0] u_y = torch.autograd.grad(u, y, grad_outputs=torch.ones_like(u), create_graph=True)[0] u_xx = torch.autograd.grad(u_x, x, grad_outputs=torch.ones_like(u_x), create_graph=True)[0] u_yy = torch.autograd.grad(u_y, y, grad_outputs=torch.ones_like(u_y), create_graph=True)[0] residual = -u_xx - u_yy - f_source(x, y) return torch.mean(residual ** 2)两次调用autograd.grad都必须写create_graph=True。第一次求u_x时开图是为了第二次能继续求u_xx;第二次求u_xx时开图是为了残差能反向传播到网络权重。任何一次漏掉,二阶导相关的那一项梯度都会静默丢失。
f_source是方程右端项,测试时常用常数 f=1,能看出解在中心处下凹。边界损失单独写一个函数:把四条边界的点集拉出来,计算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(), lr=0.1, max_iter=20) 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 是否单调下降,而不是看迭代次数。
参数设置可以直接参照这张表:
| 问题类型 | 网络结构 | 采样策略 | 优化器 |
|---|---|---|---|
| 一阶 ODE | 3 层 × 50,Tanh | 200~500 点,固定或重采样 | Adam lr=1e-3 |
| 二阶 ODE / 低维稳态 PDE | 4 层 × 80,Tanh | 1000~3000 点,重采样 | Adam 热身 + LBFGS |
| 含时间 PDE / 强非线性 | 5~6 层 × 100,Tanh 或 sin | 5000+ 点,残差自适应 | Adam 热身 + LBFGS |
含时间的 PDE 不需要特殊的网络结构,把时间 t 当普通输入特征和空间坐标拼在一起,卷积都不用加,损失函数按原方程写就行。
5. PINN 训练避坑与排查:五个让我翻车的真实原因
PINN 训练不像普通监督学习那么直观,网络翻车时你很难判断是哪一环出了问题。下面五条是按出现频率排的,每条都按现象、原因、解决的顺序写,都是我实际踩过且找到明确根源的问题。
5.1 PDE 损失不下降,只有初值项在动
现象:训练几千步,pde_loss纹丝不动,只有ic_loss在下降;或者 loss 整体下降,但解函数完全不符合方程形态。
原因:最常见的是激活函数用了 ReLU,二阶导恒为零,PDE 残差变成一个常数,梯度算不出来。另一个高频原因是torch.autograd.grad里漏了create_graph=True,一阶导以上的计算图断开,PDE 项反向传播不了。
解决:换成 Tanh 激活,并且所有autograd.grad调用都显式写create_graph=True。排查技巧是打印loss_pde.item()和loss_ic.item()两个分量,而不是只看总 loss,一眼就能看出哪项在偷懒。
5.2 训练 loss 很小,换一组点验证却对不上
现象:训练时采样点上都对,loss 压到 1e-6;把验证点取在训练点之间或区间边缘,误差突然涨到 1e-2。
原因:训练点固定且分布规则时,网络本质上在记忆这组坐标。两个训练点之间的行为只有方程残差在约束,残差对每个点贡献不均匀,网络会在采样稀疏区域“偷懒”。
解决:每个 step 重新随机采样,别用固定网格。验证时也要用全新点算残差,不要看训练 loss:
x_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,增加采样点密度;激活换回 Tanh;LBFGS 的 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_outputs=torch.ones_like(u), create_graph=True)[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 从“玄学”变成工具的唯一路径。希望帮到你。
本文还有配套的精品资源,点击获取