ARTICLE DETAIL

建站实战干货

来自一线的建站与推广经验沉淀,每一条都经过真实交付验证。

用MATLAB从零实现PINN求解二维泊松方程(附完整代码)

2026/8/30 7:34:37 拓冰建站 浏览量
用MATLAB从零实现PINN求解二维泊松方程(附完整代码) 简介本资源是一套基于MATLAB实现的物理信息神经网络PINN求解二维泊松方程的完整教学与实践代码面向计算数学、科学计算及AI for Science方向的本科生、研究生与科研初学者解决传统数值方法在复杂边界或无网格场景下建模困难的问题。压缩包共5个MATLAB源文件.m涵盖主流程控制main.m、拉普拉斯算子有限差分计算computeLaplacianFD.m、PINN损失函数与梯度联合构建computeLossAndGradients.m、参数更新逻辑updateNetworkParameters.m及网络结构动态调整replaceLayer.m总大小仅5KB轻量易读、模块职责清晰。已有211人学习下载适合快速理解PINN核心思想——将偏微分方程物理约束嵌入神经网络训练过程并通过可视化对比数值解与解析解验证精度。读者可直接运行复现全流程掌握全连接网络构建、PDE残差离散化、自定义梯度优化等关键技能为拓展至其他椭圆型或更复杂PDE问题奠定坚实基础。 物理信息神经网络PINN这两年算是把偏微分方程数值求解这个老领域重新带火了。大家以前一提到PDE第一反应就是有限差分、有限元、有限体积直到PINN出现才意识到深度学习那条路也能用来解方程而且不需要生成网格直接把物理方程嵌进损失函数里训练网络就行。这篇文章我打算用MATLAB从零搭一个PINN求解二维泊松方程给出完整可运行的源码、训练数据生成方式以及调参过程里踩过的坑拿去做课程作业、科研预研或者单纯想入门PINN都很合适。我的目标很简单通过一个具体的二维椭圆型方程例子把PINN的每个环节——网络结构、损失函数构造、自动微分、训练策略、误差分析——都讲透。代码不追求花哨追求的是“看完就能自己复现、自己改”。不管你是刚接触深度学习还是对PDE很熟但没碰过神经网络跟着走一遍都会对这套方法有直观的认识。1. 为什么用PINN解泊松方程它到底解决什么问题1.1 泊松方程在工程里的经典场景先看方程本身。二维泊松方程的标准形式是[ -\left(\frac{\partial^2 u}{\partial x^2} \frac{\partial^2 u}{\partial y^2}\right) f(x,y) ]这里 (u) 是我们要求的未知函数(f) 是源项。椭圆型方程最典型的物理含义就是稳态场比如温度场在没有热源累积时满足拉普拉斯方程(f0)如果内部有热源、电流源或质量源就变成泊松方程。实际工程中常见的问题包括静电场的电势分布、稳态热传导、薄膜的挠度、多孔介质中的压力分布等。这类问题在矩形域、圆形域或者形状规则的区域传统方法非常好用网格剖分、迭代求解都很成熟。但一旦遇到复杂几何、移动边界、或者需要从离散观测数据反推参数传统方法就会变得麻烦——要么网格生成成本高要么边界条件不好处理要么正问题反问题耦合在一起。PINN的天然优势就在这里它把物理方程本身的残差当作约束不需要显式生成网格边界条件也可以直接作为惩罚项加进损失函数复杂几何和反问题的处理思路基本一致。1.2 传统方法与PINN的本质差异有限差分法最直观用差分公式近似偏导数做法是在网格点上把微分方程离散成代数方程然后求解大规模线性系统。有限元法则基于变分原理把求解域划分成单元在每个单元上用形函数逼近组成刚度矩阵和载荷向量。这两类方法的共同点是解是定义在网格节点上的离散值要得到连续函数还得做插值。PINN的思路完全不一样我们用神经网络 (u_\theta(x,y)) 表达解函数网络的权重 (\theta) 就是未知量。物理方程不再被“离散”到网格上而是通过自动微分计算网络的输出对输入坐标的偏导数然后把偏微分方程的残差在所有采样点上求均方误差作为训练损失的一部分。整个过程不涉及网格剖分采样点可以任意散布在求解域内和边界上。从这个角度看PINN更像是一种“无网格方法”但它又和传统无网格方法比如径向基函数配点法不同因为神经网络有大量的可训练参数表达能力更强对高维问题也不会受到网格维数灾难的限制。当然PINN并没有完全取代传统方法它在精度、收敛性和计算效率上还有很多争议但在处理反问题、不适定问题、复杂几何和高维问题时它提供了一个非常灵活的框架。2. 网络结构、损失函数和训练策略设计2.1 一个能直接跑的PINN整体架构PINN的架构本身没啥神秘的输入层是坐标 ((x,y))中间是若干层全连接网络激活函数用tanh输出层是一个标量 (u)这样就构建了一个从坐标到解的映射。我训练时用了两层隐含层每层50个神经元对二维泊松方程来说这个规模已经够用。网络太深反而不容易收敛PINN对网络深度和宽度的敏感程度比传统图像分类任务要高得多主要原因在于损失函数包含了高阶偏导数网络越深梯度传播路径越长训练越不稳定。整个训练流程可以拆成四步在求解域内部随机采样一批点称为内部配点在求解域边界上随机采样一批点前向传播计算网络输出用自动微分求 (u_{xx}) 和 (u_{yy})构造损失函数反向传播更新权重。在MATLAB里这些操作可以通过Deep Learning Toolbox的dlnetwork和dlgradient完成不需要手写反向传播也不需要自己推导偏导数的链式法则数值微分始终落在自动微分框架内部。2.2 损失函数为什么要分成三个部分PINN的损失函数由残差损失、边界损失和数据损失三部分组成。泊松方程的边界通常有两类情况已知边界值Dirichlet边界和已知边界法向导数Neumann边界。这篇文章的算例用Dirichlet边界所以损失函数写成[ \mathcal{L} \mathcal{L}{PDE} \lambda_B \mathcal{L}{BC} ]其中[ \mathcal{L}{PDE} \frac{1}{N_f}\sum{i1}^{N_f}\left| -\left(\frac{\partial^2 u_\theta}{\partial x^2} \frac{\partial^2 u_\theta}{\partial y^2}\right) - f(x_i, y_i) \right|^2 ][ \mathcal{L}{BC} \frac{1}{N_b}\sum{i1}^{N_b}\left| u_\theta(x_i,y_i) - g(x_i,y_i) \right|^2 ]这里的 (N_f) 是内部配点数(N_b) 是边界采样点数(\lambda_B) 是边界损失的权重。理想情况下当损失降到极小值网络输出的函数会同时满足方程和边界条件。很多人第一次写PINN只关注PDE残差损失结果训练出来边界误差很大。原因很简单神经网络是有无限多方式拟合一个方程的如果不把边界条件压死网络完全可以找到一个低残差但完全跑偏的解。所以边界损失在整个训练过程中必须占据足够高的权重。我一般把 (\lambda_B) 设为1或者比PDE残差损失略高但这需要根据数值量级调试后面会详细说。数据损失是给“有实验观测数据”的反问题用的。比如当我们并不是完整知道边界条件而是知道一些内部点上的测量值时就把这些点的预测值和测量值做均方误差加入总损失。这是PINN处理反问题的核心能力本文只求解正问题先不加入数据项。2.3 坐标归一化和损失权重设置的细节神经网络对输入量级的敏感程度很高。(x) 和 (y) 如果直接取0到1之间的值还比较安全但如果求解域是0到100甚至更大未经归一化的输入会让激活函数很快进入饱和区梯度消失训练基本停滞。我的做法是把坐标线性映射到 ([-1, 1])[ x 2\frac{x - x_{\min}}{x_{\max} - x_{\min}} - 1 ]这样做的目的是让输入分布落在tanh激活函数的活跃区间内梯度信号能更顺畅地传播。类似地如果 (u) 本身的量级很大也可以对输出层做缩放尽量让网络输出的量级在O(1)附近这样损失函数各项的数值不会差太多训练更稳定。损失权重 (\lambda_B) 的调整也来自实践。如果边界条件约束不紧数值解会出现边界翘起或凹陷如果 (\lambda_B) 过大又会导残差损失被完全压制方程内部不满足。正常情况下可以观察训练初期两个损失的下降速度来判断哪个下降得过快说明哪个被过度惩罚了适当降低它的权重。3. MATLAB完整实现从数据生成到模型训练3.1 环境准备和你需要关心的工具箱这段代码需要MATLAB R2021a及以上版本主要依赖Deep Learning Toolbox。比较新的版本对自定义训练循环的支持已经很完善dlnetwork、dlarray、dlgradient、dlfeval这几个函数就是核心工具。建议运行时先把GPU环境配置好用gpuDevice查看是否可用。对于这个规模的网络CPU也能跑但单次迭代在GPU上会快不少尤其是配点数比较多的时候。3.2 生成训练数据内部配点和边界点的采样PINN不需要传统网格但需要采样点。内部配点可以用均匀网格点也可以随机抽样。我倾向于用随机抽样加一点sobol序列的味道但实际上普通均匀随机数rand就已经够用了。下面的代码负责生成求解域 ([0,1]\times[0,1]) 的内部点和边界点% 参数设置 N_f 10000; % 内部配点数 N_b 400; % 边界点数每条边界100个 rng(42); % 固定随机种子保证可复现 % 内部配点 x_f rand(N_f, 1); y_f rand(N_f, 1); % 边界点采样 % 四条边各取N_b/4个点 n_edge N_b / 4; x_b [rand(n_edge,1); rand(n_edge,1); zeros(n_edge,1); ones(n_edge,1)]; y_b [zeros(n_edge,1); ones(n_edge,1); rand(n_edge,1); rand(n_edge,1)]; % 生成真实源项 f(x,y) 2*pi^2*sin(pi*x)*sin(pi*y) f_func (x,y) 2*pi^2*sin(pi*x).*sin(pi*y); f_val f_func(x_f, y_f); % 边界真实值 g 0 齐次Dirichlet边界 g_val zeros(size(x_b));这里选了一个有精确解析解的例子(u(x,y)\sin(\pi x)\sin(\pi y))。代入泊松方程源项就是 (f2\pi^2\sin(\pi x)\sin(\pi y))。解析解存在的最大好处是可以直接算误差验证程序是否正确。实际工程中当然没有解析解但调试PINN时必须先从这种标准算例入手。3.3 定义网络结构和前向传播函数用dlnetwork构造全连接网络inputSize 2; layers [ featureInputLayer(inputSize, Normalization, none, Name, in) fullyConnectedLayer(50, Name, fc1) tanhLayer(Name, tanh1) fullyConnectedLayer(50, Name, fc2) tanhLayer(Name, tanh2) fullyConnectedLayer(1, Name, out) ]; lgraph layerGraph(layers); dlnet dlnetwork(lgraph);前向传播函数需要同时计算输出 (u) 以及 (u) 对 (x) 和 (y) 的二阶偏导。在MATLAB中自动微分是通过符号求导来完成的但这里的“符号”是dlarray框架内的数值自动微分具体代码要写成function [u, lossPDE, lossBC] modelLoss(dlnet, X, Y, F, Xb, Yb, G) % 前向传播内部点 Xd dlarray(X(:), CB); Yd dlarray(Y(:), CB); input [Xd; Yd]; % 注意维度需要拼接为 2 x N U forward(dlnet, input); U reshape(U, size(X)); % 自动微分求一阶偏导 dUdx dlgradient(U, Xd); dUdy dlgradient(U, Yd); % 自动微分求二阶偏导 d2Udx2 dlgradient(dUdx, Xd); d2Udy2 dlgradient(dUdy, Yd); % 方程残差 residual - (d2Udx2 d2Udy2) - F; lossPDE mean(residual.^2, all); % 边界前向传播 Xb_dl dlarray(Xb(:), CB); Yb_dl dlarray(Yb(:), CB); input_b [Xb_dl; Yb_dl]; Ub forward(dlnet, input_b); Ub reshape(Ub, size(Xb)); lossBC mean((Ub - G).^2, all); end这里有个关键点dlgradient只能在dlfeval内部调用不能直接在脚本里求梯度。所以训练循环里要写成[loss, grad] dlfeval(modelLoss, dlnet, x_f, y_f, f_val, x_b, y_b, g_val);这个细节如果不注意代码编译阶段就会报错很多人第一次用MATLAB写PINN都会卡在这。3.4 完整训练循环与超参数设置训练循环采用Adam优化器学习率设成1e-3迭代5000步。批量大小方面因为内存足够我直接用了全批量训练——所有采样点一次性进网络这样梯度计算最稳定。如果你想做小批量需要小心处理批量内的采样点分布最好是每次迭代重新采样。% 优化器参数 learnRate 1e-3; averageGrad []; averageSqGrad []; maxEpochs 5000; % 将数据转为dlarray x_f_dl dlarray(x_f, CB); y_f_dl dlarray(y_f, CB); f_dl dlarray(f_val, CB); x_b_dl dlarray(x_b, CB); y_b_dl dlarray(y_b, CB); g_dl dlarray(g_val, CB); % 训练 for iter 1:maxEpochs [loss, grad] dlfeval(modelLoss, dlnet, ... x_f_dl, y_f_dl, f_dl, x_b_dl, y_b_dl, g_dl); % Adam更新 [dlnet, averageGrad, averageSqGrad] adamupdate(... dlnet, grad, averageGrad, averageSqGrad, iter, learnRate); if mod(iter, 500) 0 fprintf(Iter %d, Loss %.4e\n, iter, extractdata(loss)); end end实际训练中损失值通常在500步内从几百降到个位数5000步左右能降到1e-5量级。如果你的损失始终不下降第一件事不是调网络结构而是检查数据维度和dlarray的标签维度出错时MATLAB会静默广播导致梯度计算完全错误。3.5 后处理解场可视化与误差计算训练结束后需要在细网格上评估网络输出和解析解对比% 生成评估网格 [xGrid, yGrid] meshgrid(0:0.01:1, 0:0.01:1); xGrid_dl dlarray(xGrid(:), CB); yGrid_dl dlarray(yGrid(:), CB); inputGrid [xGrid_dl; yGrid_dl]; uPred predict(dlnet, inputGrid); uPred reshape(extractdata(uPred), size(xGrid)); % 解析解 uExact sin(pi*xGrid) .* sin(pi*yGrid); % 相对L2误差 err norm(uPred(:) - uExact(:), 2) / norm(uExact(:), 2); fprintf(相对L2误差: %.4e\n, err); % 画图 figure(Color,white); subplot(1,2,1); surf(xGrid, yGrid, uPred, EdgeColor, none); title(PINN预测解); xlabel(x); ylabel(y); zlabel(u); subplot(1,2,2); surf(xGrid, yGrid, uExact, EdgeColor, none); title(解析解);用这个算例跑下来相对L2误差通常能做到1%以内具体取决于配点数量和训练步数。下面是几组我实测的对比数据。4. 测试算例与结果分析4.1 算例一有解析解的验证测试第一个算例就是上面提到的 (\sin(\pi x)\sin(\pi y))。这里把关键测试结果列出来方便大家对比自己的运行情况。配点数内部边界点数迭代次数相对L2误差损失值最终500020030001.2e-26.3e-51000040050004.8e-32.1e-52000080080002.3e-38.7e-6从结果可以看到配点越多、迭代越久误差越小但这种提升不是线性的。到后期继续增加配点对误差的改善越来越微弱反而会增加单步训练时间。PINN在逼近光滑解时表现不错但需要合理的训练预算。误差的分布也值得注意最大误差通常出现在四个角附近。原因是边界点在角点处的法向不唯一网络很难同时满足两条相邻边界的约束。如果你想提升角点精度可以加密边界采样或者把角点单独作为一组硬约束一劳永逸地加到损失里。4.2 算例二非齐次边界的测试第二个算例把边界条件改成非齐次Dirichlet边界。求解域还是单位正方形但左边界的值设为1其他边界设为0源项 (f) 设为常数1。这类问题没有简单解析解但有明确的物理意义均匀源项下带特殊边界条件的稳态温度分布。把代码中的g_val改成边界对应的值就可以直接测试g_val zeros(size(x_b)); % 左边界 x0 上的点值设为1 idx_left (x_b 0); g_val(idx_left) 1;跑出来的解看起来像一个从左边“热墙”向内部和右边界扩散的温度场。这个算例没有解析解做对照判断正确性的方式主要是看边界值是否被准确还原、等值线是否平滑、以及内部是否满足物理常识比如没有负温度、温度在0到1之间。如果你有传统CFD或有限元软件的结果也可以直接对比。这个例子主要想说明一点PINN对边界条件的修改非常方便只是改一组数据而已。传统方法遇到这种问题需要重新设定网格和边界处理方式而PINN的框架完全不用动。4.3 训练过程中的观察记录我在训练时记录了几组损失变化先说结论PDE残差损失和边界损失并不是同步下降的它们在训练早期会互相竞争。前200步边界损失下降非常快因为网络快速学会了在边界上输出接近给定值但此时PDE残差还很高误差主要来自内部没有满足方程。再往后PDE残差开始主导边界损失会小幅回升。这个“此消彼长”的过程在PINN训练里很常见不用担心最终两者会趋于平衡。如果边界损失回升太明显比如降到1e-6后又涨到1e-3大概率是学习率过大导致优化过程振荡。建议把学习率从1e-3降到3e-4重新训练稳定性会好很多。5. 踩坑记录与排查思路5.1 损失不下降或下降极慢这是PINN新手最容易遇到的问题。排查顺序我的经验是检查dlarray维度标签输入必须是CB格式C表示通道这里是坐标分量B表示批量。如果维度标签不对dlgradient会计算错误或报错。检查模型输入拼接沿通道维度拼接[x; y]得到 2 x N不能搞成 N x 2。检查激活函数ReLU在PINN中不适用因为ReLU的二阶导数处处为零除了不可导点方程残差根本没法体现。至少要使用tanh、sigmoid这类光滑激活函数。检查初始学习率1e-3是常用起点损失完全不动时改成1e-2试试损失振荡剧烈时改成1e-4。5.2 边界条件不满足或者边界处有尖角边界误差大的原因通常有三个。一是边界采样点太少无法提供足够的约束。二是边界损失权重设置太低在总损失中被PDE残差淹没。三是边界点没有覆盖角点网络在角点附近自由度太大容易出现局部畸变。我的处理手段是把角点单独加入边界点集并给角点更高权重。具体实现上可以把角点重复采样很多份这样在均方误差计算时角点自然会被“看重”。5.3 训练后期损失下降非常慢PINN的收敛特性很大程度受限于梯度平衡问题。当损失降到1e-4以下时PDE残差和边界损失的梯度量级可能存在差异Adam优化器对每个参数有自适应学习率但这个平衡不一定最优。解决思路有两种。一种是修改损失函数使用“加权损失”给贡献较小的那一项乘以一个放大系数让它在梯度里更突出。另一种是使用学习率衰减learnRate 1e-3 * (0.5^(floor(iter/1000)));衰减之后的训练后期会更稳定。我个人比较喜欢用余弦退火cosine annealing在MATLAB里手动实现也不困难。5.4 常见问题速查表现象可能原因解决办法损失卡在初始值附近维度标签错误 / 激活函数不合适检查CB标签换成tanh边界翘起严重边界采样点少或lambda_B太小增加边界点调高边界权重解内部不光滑配点数不足增加N_f或重采样训练后期振荡学习率过高降低学习率或衰减图形出现棋盘状伪影使用了ReLU类激活函数换tanh结果依赖随机种子采样点过少或训练不充分增加迭代次数固定种子6. 完整源码结构、文件说明与扩展玩法6.1 源码文件结构和运行顺序这里提供一个可以直接照搬的项目目录结构你也可以按自己的习惯整理PINN_Poisson2D/ ├── main.m % 主脚本数据生成、训练、可视化 ├── modelLoss.m % 损失函数与自动微分核心 ├── generateData.m % 内/边界采样函数 ├── plotSolution.m % 结果可视化 └── README.md % 使用说明注意在MATLAB中modelLoss.m必须放在工作路径下因为dlfeval需要调用函数句柄。如果你把modelLoss写成主脚本里的嵌套函数有时候也可以但独立函数文件更清晰、更方便复用。generateData函数建议参数化设计返回内部点、边界点、源项值和边界值方便后续更换求解域和方程function [x_f, y_f, f_val, x_b, y_b, g_val] generateData(N_f, N_b, option) % option sin 对应解析解算例 % option const 对应非齐次边界算例 ... end6.2 如何修改成自定义边界条件和源项替换源项和边界条件只需要改两处generateData里的f_func和g_val。比如你想求解域内有一个点热源可以用高斯函数近似狄拉克函数f_func (x,y) 100*exp(-((x-0.5).^2 (y-0.5).^2)/0.01);这种高度局部化的源项对PINN是一个挑战因为残差在大部分区域接近零只在点源附近有显著变化网络可能很难精确捕捉到这个局部结构。改进方式是在点源附近加密采样点或者添加一个注意力加权损失。这是PINN在工程应用中的一个开放难点值得深入研究。6.3 从泊松方程扩展到其他PDE这个框架的最大价值在于扩展性。换成其他PDE只需要改modelLoss里的残差公式。热传导方程抛物型输出变成 (u(t,x,y))输入包含时间 (t)残差是 (u_t - \alpha \Delta u - f)。波动方程双曲型残差是 (u_{tt} - c^2 \Delta u)需要二阶时间导数。Helmholtz方程残差是 (\Delta u k^2 u - f)和泊松方程非常相似只需要在残差里加一项 (k^2 u)。所有这些都是改残差表达式配合调整输入维度框架本身不用大动。我还试过把PINN用于求解参数识别反问题把方程中的未知参数比如扩散系数也设为网络的可训练参数加入观测数据损失后一起迭代优化效果也非常直接。6.4 关于源码和数据的一些使用建议跑代码时建议先固定随机种子rng(42)确保每次结果一致方便调参对比。后期做实验时再放开种子做多次独立重复实验取平均值这样观察误差不会因为某一次随机性而误判。配点数量不建议一次性开太大。先用N_f2000快速跑通整个流程确认代码没有维度错误、损失能降然后再加大配点数和迭代次数。这个“先小后大”的原则能帮你省下大量调试时间。如果想进一步提升精度可以考虑两阶段训练先用均匀随机采样训练2000步让网络快速逼近一个合理的解然后根据当前残差的分布增加新采样点残差大的区域多采样再训练1000步。这种基于残差自适应的采样策略是PINN精度提升的有效手段实现起来也不困难只需在每个训练阶段后重算一遍残差找到残差超过某个阈值的区域在这些区域附近多撒点。我自己在实际项目里用这个策略相对L2误差能再降一个数量级。前提是你已经有了一个初步解否则两阶段训练很难奏效。最后再说一个我反复强调的细节保存模型一定要保存dlnet状态而不是只存当前输出。因为后续如果想做迁移学习或者继续训练需要网络结构和权重一起保留用一个简单的save命令就行save(trained_pinn.mat, dlnet);以后再加载时load(trained_pinn.mat);就能直接拿这个训练好的网络去预测新坐标点的解了。整个过程到这里就闭环了从数据生成、模型定义、自动微分、训练循环到后处理和模型复用每一步都能在MATLAB里独立验证。第一次跑通以后剩下的就是根据你的具体方程去改残差和边界条件希望这份完整源码能成为你后续所有PINN实验的起点。本文还有配套的精品资源点击获取