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 代码包

PINN（Physics-Informed Neural Networks）求解微分方程的门槛，从来不在理论推导，而在把方程写进 Python 代码的那一刻。这套资源把十个以上 PINN 场景做成了直接可跑的 Notebook：从最简单的一阶 ODEy'=sin(π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后缀的文件是作者调试时留下的备份，可以忽略，但别删除——万一改了主文件跑崩了，备份就是后悔药。

先看用例层，按从易到难整理了一份清单：

看到Δu=2, x∈[-1,1]这个文件名别奇怪，它其实是把一维 Poisson 方程u''=2配了一个 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 接口，太新的版本在某些老电脑上会出现编译指令集的警告，不影响跑通但输出很吵。装完依赖后启动 Jupyter：

jupyter notebook

浏览器打开 Notebook 列表后，第一遍建议先跑y'=sin(πx)cos(πx), x∈[-1,1], y(-1)=0.ipynb。这个文件是整套资源里最简的闭环：一阶常微分方程、解析解可以手写验证、网络只需要两三层。先把最小闭环跑通，再进 Poisson 和扩散方程，排查问题的范围会小很多。

3.2 以 Poisson 方程为例：损失函数与训练主循环

Poisson 方程通常写作-∇²u = f。为了和文件名里的Δu = 2直接对应，下面代码改用u'' = f的形式，即在Δu=2这个例子里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(persistent=True) 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必须设置persistent=True，因为后面如果还要算u_yy就需要对同一个梯度带调用多次 gradient。残差是u'' - f，当网络输出收敛到真解时，残差处处趋近于零，均方误差也趋近于零。

参数说明：x_f是内部配点，通常用均匀或随机采样生成，数量取 1000 到 5000 之间；f由具体问题决定，Δu=2例子里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_rate=1e-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(f"step {step}, loss {loss.numpy():.2e}")

跑完 Adam 阶段再切 L-BFGS：PyTorch 里直接用torch.optim.LBFGS，需要把整个 loss 计算逻辑包进 closure；TensorFlow 里一般借助 scipy 的scipy.optimize.minimize(method='L-BFGS-B')或 TensorFlow Probability 的lbfgs_minimize接口，把trainable_variables拍平成向量传入。切过去之后如果损失还在缓慢下降，说明方向是对的；如果损失反而上升，说明 Adam 阶段跑太长已经过拟合了配点，这时候要减少 Adam 步数或者增加配点数量。

我自己跑这套资源有一个固定动作：每换一个新方程，先把配点数降到很小的值（比如 100 个点）快速跑通流程确认没报错，再恢复完整配点数去调精度。这个小习惯能省掉至少一半的排错时间。

4. 边界条件与特殊域：五种写法、圆盘极坐标与方程组

4.1 五类边界条件的损失函数设计

读这套代码的 Notebook 清单时你会发现，光 Poisson 方程就有五个变体：Dirichlet、Dirichlet_Robin、Dirichlet_Periodic、Dirichlet_Neumann、以及纯 Dirichlet。作者其实是在用同一套方程演示边界条件的五种写法，这是这份资源里最值得反复对比的部分。

先把五类边界条件的定义和损失写法放一张表：

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(persistent=True) 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 降得很漂亮，相图却完全收敛到错误吸引子上。从那以后，我每次接触一个新方程，都会先拿一个带解析解的特例把残差、边界、优化器全部跑通，再换真实参数，这个习惯帮我避开了至少一半的调参时间。希望帮到你。

本文还有配套的精品资源，点击获取