PINN求解微分方程：一套可直接复现的Python代码包

简介面向物理信息神经网络PINN学习者与科研人员这份压缩包提供了一套完整的Python实现案例覆盖常微分方程、偏微分方程以及Lorenz系统等典型问题并包含DeepXDE框架的泊松方程示例帮助读者理解将物理定律嵌入神经网络损失函数的核心思想。资源共26个文件以17个ipynb案例笔记为主配套3个py脚本用于模型与几何定义另有4个zbak备份文件、1个png结果图及1个md说明文档压缩包仅891KB轻量易用。已有217人学习适合希望快速上手PINN求解流程的初学者和需要参考代码框架的进阶用户。通过案例可掌握数据准备、网络结构设计、损失函数构建、Adam或LBFGS优化训练与结果分析等完整步骤并对边界条件、初值条件及不同方程类型的建模差异形成直观认识可供教学与实战参考。1. PINN 求解微分方程一套能直接复现的 Python 代码包PINNPhysics-Informed Neural Networks求解微分方程的门槛从来不在理论推导而在把方程写进 Python 代码的那一刻。这套资源把十个以上 PINN 场景做成了直接可跑的 Notebook从最简单的一阶 ODEysin(πx)cos(πx)到带 Dirichlet、Neumann、Robin、周期混合边界的 Poisson 方程再到圆盘域上的 Laplace 方程、Lorenz 混沌系统甚至包含需要四阶导数的 Euler Beam 算例。如果你卡在理论懂了但代码跑不通或自己的 PINN 不收敛这份代码包可以直接抄作业——方程、边界条件、采样区间都写在文件名里照着改就能迁移。适合做课题验证的研究生、做仿真预研的工程师也适合把 PINN 当基线方法的算法同学。2. 拆开压缩包从 PDE.py 到 model.py 的模块分工2.1 文件地图哪些是演示用例哪些是核心库先别急着双击 Notebook花五分钟搞清楚目录结构能省下后面大量排错时间。这份压缩包里的文件分三类以方程命名的 .ipynb 是演示用例每个文件对应一个完整的 PINN 求解任务model.py、PDE.py、geometry.py三个 .py 文件是核心实现被所有 Notebook 复用.zbak后缀的文件是作者调试时留下的备份可以忽略但别删除——万一改了主文件跑崩了备份就是后悔药。先看用例层按从易到难整理了一份清单文件对应的数学问题学习价值ysin(πx)cos(πx), x∈[-1,1], y(-1)0.ipynb一阶 ODE 边值问题理解 PINN 最小闭环y-y-y2x, x∈[-5,5], y(-5)1, y(5)5.ipynb二阶 ODE 边值问题二阶导自动微分写法Δu2, x∈[-1,1], u(-1)0, u(1)4.ipynb一维 Poisson Neumann边界损失里加导数Poisson_equation_Dirichlet.ipynb系列二维 Poisson 多种边界组合边界条件设计核心场景Laplace_equation_on_a_disk.ipynb圆盘域 Laplace极坐标域采样Diffusion equation.ipynb扩散方程 u_t αu_xx含时间的 PDEA_simple_ODE_system.ipynb常微分方程组多输出网络Lorenz system.ipynbLorenz 混沌系统方程组的经典压力测试Euler Beam.ipynb四阶梁方程高阶导训练避坑PDAE_system.ipynb偏微分代数方程约束耦合场景DeepXDE/用第三方库求解 Poisson与自写代码交叉验证Method Testing/Jacobian-Hessian 两种求导法二阶导优化技巧看到Δu2, x∈[-1,1]这个文件名别奇怪它其实是把一维 Poisson 方程u2配了一个 Neumann 边界u(1)4和一个 Dirichlet 边界u(-1)0这个组合在物理上对应一端固定、一端给边界流量的杆件热传导问题。作者把方程直接写在文件名里是为了方便检索但也意味着文件名里带了很多数学符号在 Windows 上解压偶尔会出现乱码。根目录下还有一张robin plot.png是作者跑 Robin 边界算例后导出的预测与解析解对比图可以作为你复现结果的参照。2.2 三个核心模块model.py、geometry.py、PDE.py 各管什么用例层之下是三个被复用的核心模块职责划分非常清晰这也是这套代码最值得抄的部分。model.py负责神经网络本身。典型实现是输入坐标点x一维问题输入是[N,1]二维问题是[N,2]经过若干层全连接层和激活函数后输出预测解u。PINN 里激活函数几乎默认用tanh原因很简单它二阶以上导数连续且非零而 ReLU 的二阶导恒为零根本无法用于 Poisson 这类二阶方程。网络深度一般 4 到 8 层宽度 20 到 100 之间具体看问题复杂度——Lorenz system.ipynb这类混沌系统往往需要更宽的网络而单纯的一阶 ODE 三层就够。geometry.py负责计算域。它生成两类东西内部配点collocation points和边界点。内部配点用来计算方程残差边界点用来计算边界条件损失。采样方式有均匀网格和随机采样两种PINN 里通常用拉丁超立方或简单随机采样配点数量级在一千到一万之间。这里有个容易被忽视的细节边界点必须被网络看见足够多次否则边界条件就是软约束网络很容易找到一个满足内部方程但在边界上漂移的解。PDE.py负责把方程翻译成损失函数。它把网络输出、模型对坐标的导数、方程右端项组合成残差表达式再对残差取均方误差。三个模块的调用关系是Notebook 里先实例化 geometry 拿到采样点再实例化 model 拿到网络最后把两者交给 PDE 构建损失函数。这个解耦设计的好处是换方程只需要改 PDE.py 里的残差表达式换几何域只需要改 geometry.py网络结构基本不用动。2.3 PINN 的完整求解链路把三个模块串起来一次标准的 PINN 训练循环是这样走的内部采样点x_f和边界采样点x_b输入网络前向传播得到预测值u(x_f)和u(x_b)对u关于输入坐标求一阶导和二阶导代入 PDE 残差表达式得到残差损失项边界条件同样写成损失项比如 Dirichlet 边界就是(u(x_b) - g(x_b))²总损失是残差项、边界项、初值项时间问题里的加权和然后对网络参数做梯度下降。数据域这一层作者单独放了一个data_domain.ipynb它专门演示采样点如何生成、如何可视化。实际跑其他 Notebook 之前先把这个过一遍能直观理解域和边界在 PINN 里到底长什么样——很多时候你觉得自己在调网络其实是在调采样点分布。这套结构的最大价值在于所有 Notebook 共享同一套底层代码因此可以把每个算例看作一个参数指纹。当你把某个 Notebook 跑通后换方程时只需要复制一份改动几何域、残差表达式和边界损失三处剩下的训练流程完全复用。下一章就按这个思路从环境准备开始完整带一遍。3. 动手复现从一阶 ODE 到 Poisson 方程的完整流程3.1 环境准备Python 解释器与依赖版本这份资源下载解压后依赖面非常窄属于开箱即用的类型。核心需要 numpy数组运算与采样、matplotlib绘图、以及一个带自动微分的深度学习框架。Notebook 里用到的自动微分能力TensorFlow 2.x 的 GradientTape 和 PyTorch 的 autograd 都能覆盖我以 TensorFlow 2.x 的写法为例说明PyTorch 用户只需要把求导接口替换成torch.autograd.grad网络定义换成nn.Module代码骨架不变。安装依赖用 pip 一条命令解决pip install numpy matplotlib tensorflow jupyter版本建议Python 用 3.9 到 3.11 之间TensorFlow 用 2.10 到 2.15 之间。太老的 TensorFlow 1.x 没有 GradientTape 接口太新的版本在某些老电脑上会出现编译指令集的警告不影响跑通但输出很吵。装完依赖后启动 Jupyterjupyter notebook浏览器打开 Notebook 列表后第一遍建议先跑ysin(πx)cos(πx), x∈[-1,1], y(-1)0.ipynb。这个文件是整套资源里最简的闭环一阶常微分方程、解析解可以手写验证、网络只需要两三层。先把最小闭环跑通再进 Poisson 和扩散方程排查问题的范围会小很多。3.2 以 Poisson 方程为例损失函数与训练主循环Poisson 方程通常写作-∇²u f。为了和文件名里的Δu 2直接对应下面代码改用u f的形式即在Δu2这个例子里f恒等于 2。Poisson_equation_Dirichlet.ipynb把 Dirichlet 边界条件单独做成一个损失函数和残差损失组合起来就是完整的 PINN 目标函数。下面这段代码是 PINN 损失函数的核心实现和压缩包内PDE.py的组织方式一致import numpy as np import tensorflow as tf def pde_loss(net, x_f, f): # x_f: 内部配点, shape [N, 1] with tf.GradientTape(persistentTrue) as tape2: tape2.watch(x_f) with tf.GradientTape() as tape1: tape1.watch(x_f) u net(x_f) # 网络前向传播 u_x tape1.gradient(u, x_f) # 一阶导 u_xx tape2.gradient(u_x, x_f) # 二阶导 residual u_xx - f(x_f) # 方程残差: u - f 0 return tf.reduce_mean(tf.square(residual))代码逻辑说明net是上一章 model.py 里定义的全连接网络f是方程右端项函数。第一层 GradientTape 用来求一阶导u_x第二层在u_x的基础上再求一次导得到u_xx这就是二阶导数的标准写法。tape2必须设置persistentTrue因为后面如果还要算u_yy就需要对同一个梯度带调用多次 gradient。残差是u - f当网络输出收敛到真解时残差处处趋近于零均方误差也趋近于零。参数说明x_f是内部配点通常用均匀或随机采样生成数量取 1000 到 5000 之间f由具体问题决定Δu2例子里f恒等于 2。换方程时只需要改 residual 这一行比如换成扩散方程就写u_t - alpha * u_xx换成梁方程就写u_xxxx - q。边界条件的损失函数单独写def boundary_loss(net, x_b, u_b): # x_b: 边界采样点, u_b: 边界给定的解值 pred net(x_b) return tf.reduce_mean(tf.square(pred - u_b))把两个损失加权组合成总损失进入训练循环total_loss pde_loss(net, x_f, f) lam_b * boundary_loss(net, x_b, u_b)lam_b是边界损失权重这是 PINN 调参里最重要的超参数之一取值在 1 到 100 之间第五章会专门展开它引起的坑。3.3 训练策略Adam 起步LBFGS 收尾这套代码包里大部分 Notebook 的训练循环遵循同一个套路先用 Adam 优化器跑几千步再用 L-BFGS 精修。原因是 PINN 的损失面有大量平坦区域Adam 这类一阶方法前期收敛快但后期在平坦区挪不动L-BFGS 作为拟牛顿法利用二阶曲率信息能让损失在最后阶段再降一到两个数量级。optimizer tf.keras.optimizers.Adam(learning_rate1e-3) tf.function def train_step(): with tf.GradientTape() as tape: loss total_loss() grads tape.gradient(loss, net.trainable_variables) optimizer.apply_gradients(zip(grads, net.trainable_variables)) return loss for step in range(2000): loss train_step() if step % 500 0: print(fstep {step}, loss {loss.numpy():.2e})跑完 Adam 阶段再切 L-BFGSPyTorch 里直接用torch.optim.LBFGS需要把整个 loss 计算逻辑包进 closureTensorFlow 里一般借助 scipy 的scipy.optimize.minimize(methodL-BFGS-B)或 TensorFlow Probability 的lbfgs_minimize接口把trainable_variables拍平成向量传入。切过去之后如果损失还在缓慢下降说明方向是对的如果损失反而上升说明 Adam 阶段跑太长已经过拟合了配点这时候要减少 Adam 步数或者增加配点数量。我自己跑这套资源有一个固定动作每换一个新方程先把配点数降到很小的值比如 100 个点快速跑通流程确认没报错再恢复完整配点数去调精度。这个小习惯能省掉至少一半的排错时间。4. 边界条件与特殊域五种写法、圆盘极坐标与方程组4.1 五类边界条件的损失函数设计读这套代码的 Notebook 清单时你会发现光 Poisson 方程就有五个变体Dirichlet、Dirichlet_Robin、Dirichlet_Periodic、Dirichlet_Neumann、以及纯 Dirichlet。作者其实是在用同一套方程演示边界条件的五种写法这是这份资源里最值得反复对比的部分。先把五类边界条件的定义和损失写法放一张表边界条件数学形式损失函数写法Dirichletu(x) g(x)mean((u_pred - g)²)Neumannu(x) h(x)mean((u_x_pred - h)²)Robina·u(x) b·u(x) cmean((a·u_pred b·u_x_pred - c)²)Periodicu(x_0) u(x_1)mean((u_pred(x_0) - u_pred(x_1))²)混合组合上述任意组合各项加权求和Neumann 边界和 Dirichlet 边界的区别只在损失里用的是网络输出本身还是一阶导数值。写 Neumann 边界损失时注意必须复用求导那层的梯度带如果单独为边界点再开一个梯度带就要用同一个模型结构否则容易造成重复计算或者梯度断链。Periodic 边界是这套资源里比较容易被误解的一个。它约束的不是某个点的值等于给定值而是域的两个端点值相等。在Poisson_equation_Dirichlet_Periodic.ipynb里做法是分别对左边界和右边界采样让网络在两个边界上的输出差值的平方趋近于零。这是 PINN 实现周期条件的标准做法比强行改造网络输入特征要简单得多。Robin 边界出现在Poisson_equation_Dirichlet_Robin.ipynb。它是 Dirichlet 和 Neumann 的线性组合在传热问题里对应第三类边界条件环境对流换热。写损失的时候把线性组合表达式整体放进残差项即可不需要拆分。robin plot.png那张对比图就是这个算例的输出横轴是空间坐标两条曲线分别是解析解和 PINN 预测Robin 边界所在的那一端曲线斜率是给定的而不是函数值给定看图时能明显感觉到两种边界的约束差异。4.2 圆盘域与极坐标Laplace_equation_on_a_disk 的采样思路Laplace_equation_on_a_disk.ipynb是这份资源里几何上最有意思的一个。圆盘域如果直接用直角坐标采样边界很难精确描述所以常见做法是切换到极坐标令x r·cosθ、y r·sinθ在(r, θ)空间里圆盘变成一个规则的矩形域r∈[0, R], θ∈[0, 2π]。换坐标之后Laplace 算子要跟着换。平面极坐标下的拉普拉斯算子是Δu u_rr (1/r)·u_r (1/r²)·u_θθ这比直角坐标下的u_xx u_yy多出两个系数项而且1/r在圆心处是奇异的。实际处理时圆心点要么单独避开要么在采样时保证r最小值不为零。训练损失里对预测u求关于r和θ的导数再用上面这个表达式组合成残差几何问题就转化成了坐标变换问题。这种处理让边界采样变得极其干净——圆盘边界就是矩形域里r R的一条边直接采一条线就行。这种坐标变换算子替换的思路对不规则几何很有价值以后遇到圆环、扇形、球壳这类域都可以先把几何映射到规则域再构造损失比在原坐标系里硬算要稳健得多。4.3 从单方程到方程组ODE system 与 Lorenz 混沌系统A_simple_ODE_system.ipynb和Lorenz system.ipynb把问题从单个方程推到了方程组。方程组和单方程的区别在于输出维度网络不再输出一个标量u而是输出一个向量[u1, u2, u3]每个分量对应一个未知函数。Lorenz 系统是三个一阶常微分方程dx/dt σ(y-x)dy/dt x(ρ-z) - ydz/dt xy - βz典型的参数取σ10, β8/3, ρ28这个取法下系统进入混沌状态对初值极其敏感。用 PINN 解它的难点不在损失函数——三个方程的残差分别计算再相加即可——而在训练稳定性混沌系统一点点数值误差都会被时间演化放大所以配点的时间范围往往不能取得太大且学习率要调得比普通问题更小。作者把 Lorenz 单独放一个 Notebook本质上是拿最敏感的场景做压力测试你在自己的问题上如果遇到训练不稳可以参考它的处理方式。PDAE_system.ipynb里涉及的偏微分代数方程更复杂它同时包含偏微分方程和纯代数约束损失函数里必须显式加入代数约束项。这一类题目在传统数值方法里处理起来很绕但 PINN 天然能处理——约束本质上就是额外一项损失。这也是 PINN 相对传统数值方法真正的优势场景之一。5. 避坑记录PINN 训练不收敛时的四个高频原因5.1 损失下降的假象与收敛停滞两条典型的训练期问题踩坑记录一损失已经降到 1e-5但预测解完全不对。现象是训练日志里总损失一路下降甚至降到 1e-5 以下但把网络预测画出来一看和真实解差得离谱有时是接近零的平坦直线。原因是残差项和边界项之间的量级失衡。PINN 的总损失是各项加权和当方程右端项数值很大时比如 Poisson 方程右端项是 100 这个量级网络只要输出一个满足边界条件的常数残差项就接近零了总损失也可以很小但这个常数解显然不是真解。损失下降只说明梯度在降低不说明解在逼近真解。解决方法是先把方程做归一化处理让解和右端项的量级都落在 O(1) 附近然后分别打印残差损失和边界损失的数值如果两者差两个数量级以上就给边界损失乘一个权重系数从 10 开始尝试直到两者落在同一数量级。我一般会在训练脚本里加一段监控代码每 100 步同时打印总损失、残差损失、边界损失三个数谁失衡一目了然。踩坑记录二Adam 跑了两千步损失卡在同一个值不动。现象是训练前几百步 loss 下降很快之后进入缓慢爬行模式无论怎么调学习率损失都停在同一水平线。原因是 PINN 的损失函数是非凸的Adam 这类一阶优化器在平坦区域只有一个方向的梯度信息步长被压缩得极小看起来就是不动了。这不是代码 bug而是优化器特性。这套资源里多个 Notebook 能在最终结果上达到很高精度靠的正是 Adam 之后又切了 L-BFGS。解决方法是让 Adam 只负责把网络从随机初始化拉到一个合理的盆地跑 2000 到 5000 步后切换 L-BFGS 做精修。切换时会新构造优化器Adam 的动量状态自然就被丢弃了不需要额外处理如果切过去后损失反而上升说明 Adam 阶段过拟合了配点减少 Adam 步数或者增加配点数量即可。5.2 高阶导崩溃与备份文件干扰两条结构性坑踩坑记录三Euler Beam 这类四阶导算例梯度爆炸或 NaN。现象是跑Euler Beam.ipynb时前几百步 loss 正常某一步突然变成 NaN之后再也救不回来重启训练在同一个位置附近又会崩一次。原因是四阶导数对网络输出的细微抖动极其敏感全连接网络在训练初期权重较大时输出的高阶导数值可以轻易达到成百上千的量级平方后直接溢出ReLU 系列激活函数在这种场景下更是碰都不能碰它的二阶以上导数不是零就是不存在。解决方法是第一个办法把激活函数统一换成tanh它的所有阶导数连续且有界第二个办法把网络做宽而不是做深深度超过八层的 PINN 在训练初期的数值行为很难控制宽度从 32 加到 128 通常更有效第三个办法是降阶把四阶方程拆成两个二阶方程组用两个网络输出分别表示w和w让每一层的导数阶数降下来。这三个办法按顺序试能解决绝大多数高阶导崩溃问题。踩坑记录四改了 PDE.py 但结果没变化或者打开的是备份文件。现象是修改了PDE.py里的残差表达式重新运行 Notebook结果和修改前一模一样或者代码报错提示找不到某个函数但明明在文件里看见过。原因是.zbak文件在这份资源里干扰性很强作者在调试时把一份PDE.py备份成PDE.py.zbak部分编辑器和文件管理器默认按扩展名排序.py.zbak和.py排在一起肉眼很容易看错。另一个常见情况是 Jupyter 在启动时加载了旧的.py模块缓存修改文件后没有重启内核import 到的还是旧版本。解决方法是训练前先在终端里用ls -la确认要修改的文件真实存在且拼写正确修改完.py文件后在 Notebook 里执行import importlib; importlib.reload(PDE)强制重载或者直接重启内核。还有一个更省事的习惯把所有.zbak文件挪到一个backup/子目录里从根源上消除干扰。6. 进阶技巧Jacobian-Hessian 求导与损失权重的实操经验6.1 Jacobian-Hessian 写法一次求导拿到所有二阶项二维 Poisson 的损失函数需要u_xx u_yy。新手很容易写成对 x 求一层导数、对 y 单独再开一个梯度带这样数据在两层梯度带之间反复拷贝计算图冗余。更常见的做法是用一个 persistent 梯度带把一阶分量一次取全再分别求二阶项import tensorflow as tf def laplace_loss(net, xy): # xy: 内部配点, shape [N, 2]第一列 x第二列 y with tf.GradientTape(persistentTrue) as tape2: tape2.watch(xy) with tf.GradientTape() as tape1: tape1.watch(xy) u net(xy) # 网络输出, shape [N, 1] grads tape1.gradient(u, xy) # 一阶导, shape [N, 2] u_x, u_y grads[:, 0], grads[:, 1] u_xx tape2.gradient(u_x, xy)[:, 0] # 关于 x 的二阶导 u_yy tape2.gradient(u_y, xy)[:, 1] # 关于 y 的二阶导 return tf.reduce_mean(tf.square(u_xx u_yy))代码逻辑说明外层 tape2 里同时记录了u_x和u_y的计算路径所以能对它们再求一次导[:, 0]和[:, 1]分别取出关于 x 和 y 的二阶分量相加得到拉普拉斯算子。这样整个计算图只开两层梯度带和分别求两次一阶导相比省掉一半的自动微分开销。Method Testing 目录里的Jacobian-Hessian methods for Laplace equation.ipynb演示的正是这个对比——小配点数下两者没差别配点数上万时合并写法的训练速度优势非常明显。6.2 损失权重、配点与物理校验的收尾习惯损失权重没有万能解但有一条快速的起点路径先把所有损失项初始权重都设为 1跑 100 步看各项量级然后按让每一项对总损失的贡献在同一数量级的原则调整。残差项数值偏大就调小它的权重边界项偏小就调大一般只动边界权重就能解决绝大多数问题权重范围在 0.01 到 100 之间尝试就够不需要追求精细。配点数量和分布同样影响结果。内部配点决定残差约束的密度边界配点决定边界条件被满足的程度。常见做法是让边界点数量是内部点的十分之一到三分之一而不是各一半——边界附近解的剧烈变化需要足够多的点来约束但边界点不提供方程的物理信息堆太多对内部精度没有帮助。最后也是最容易被忽略的一步用物理量验证解而不只是看 loss。扩散方程算例里检查总质量是否守恒Lorenz 系统看相轨迹是否落在奇怪吸引子上Poisson 问题用解析解做逐点最大误差。这套资源里每一个 Notebook 都配了解析解或参考解这就是给你对比用的。我从这套代码里最深的教训是第一次跑 Lorenz 时loss 降得很漂亮相图却完全收敛到错误吸引子上。从那以后我每次接触一个新方程都会先拿一个带解析解的特例把残差、边界、优化器全部跑通再换真实参数这个习惯帮我避开了至少一半的调参时间。希望帮到你。本文还有配套的精品资源点击获取