ARTICLE DETAIL

资讯详情

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

MATLAB实现PINN求解二维泊松方程:完整源码与实战指南

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解、有限差分辅助微分算子计算等核心机制是掌握PINN原理与MATLAB工程实现的优质入门范例。 我先把话说在前头如果你只是想快速算一个二维泊松方程传统有限差分法五分钟就能写完真没必要上PINN。但这篇文章我针对的是另一类人——你已经厌倦了反复剖网格、调离散格式、换一个边界条件就重写一遍求解器的人。PINN这个思路把偏微分方程本身变成损失函数用神经网络直接拟合解函数MATLAB里从零搭一遍完整流程你会对“求解PDE”这件事有完全不同的理解。这篇文章会给你一份能直接跑通的MATLAB源码包括训练数据的生成策略、神经网络结构选择、损失函数写法、自动微分处理二阶导数的细节以及训练完之后的误差验证方法。适合有MATLAB基础、对神经网络有概念但没亲手写过PINN的读者也适合已经用Python跑过PINN、想对比看看MATLAB实现差异的人。1. 为什么泊松方程值得用PINN来解1.1 传统网格类解法的繁琐与PINN的“零网格”思路二维泊松方程是工程里最常见的椭圆型方程之一。电势分布、稳态热传导、不可压缩流体的压力场最后都会落到这个形式上。传统做法是先把求解区域剖成网格然后在每个节点上构造代数方程最后解一个大型稀疏线性方程组。思路本身不复杂但真正做起来前期网格剖分和后处理取点这两件事就够折腾的。区域稍微不规则一点网格质量的控制就成了玄学。PINN的思路完全反过来了。它不剖网格而是把解函数表示成一个神经网络输入是坐标(x, y)输出是解u的值。训练这个网络的时候损失函数由两部分构成一部分是方程残差即在采样点上计算u的偏导数并带入方程左边看它和源项差多少另一部分是边界条件残差即网络输出和已知边界的偏差。训练让这两部分同时趋近于零网络学到的函数就是方程的近似解。这个思路最吸引人的地方在于一旦框架搭好换边界条件、换源项甚至换求解区域都只是改采样点和损失函数的问题不需要重新设计离散格式。这也是为什么我最初决定在MATLAB里把完整流程跑通的原因——验证一下这个“无网格求解”在工作流上到底能省多少事。1.2 泊松方程的物理背景与本文的基准问题为了验证代码的正确性我用一个带解析解的基准问题来测试这是调试PINN最稳妥的做法。求解区域取单位正方形方程形式是-u_xx - u_yy f(x, y), (x, y) ∈ [0,1]×[0,1]边界条件取Dirichlet零边界即边界上u 0。我选的解析解是u(x, y) sin(πx) * sin(πy)把解析解代入方程左边求导可以得到对应的源项f(x, y) 2π² * sin(πx) * sin(πy)选择这个基准问题有三个原因。第一边界条件简单网络训练时边界残差部分实现容易第二解析解光滑神经网络拟合起来比较轻松适合验证框架的可行性第三正弦函数族的偏导计算简单方便我核对MATLAB自动微分算出来的二阶导数是否正确。这个基准问题对标的是物理学里的简单场景一个方形区域四边接地、内部有均匀分布电荷时的电势分布。实际工程问题往往更复杂但先把简单的跑通、把每一行代码都吃透再去加工程细节才有基础。2. MATLAB里搭PINN的三个技术支点2.1 自动微分dlarray、dlgradient、dlfevalMATLAB做PINN核心工具是Deep Learning Toolbox里的自动微分能力。三个关键函数要搞明白dlarray用来封装带计算图的数据dlgradient用来求偏导dlfeval用来在计算图环境下执行函数。PINN的损失函数里需要算u对x和y的二阶偏导。神经网络的前向传播本身是一个复杂的复合函数手动推导二阶导完全不现实。自动微分会在前向传播时记录每一步运算关系然后通过链式法则自动求出各阶导数。在MATLAB里这个过程被抽象成了dlgradient的两次调用第一次对前向传播结果求一阶导第二次对一阶导再次求导得到二阶导。一个容易踩的坑是dlgradient不能直接在普通函数里调用必须在dlfeval的上下文里执行。我第一次写的时候在普通脚本里直接调用dlgradientMATLAB直接报错提示只能在dlfeval内使用。这个限制是计算图机制决定的写代码时要把损失函数的计算封装在一个独立函数里然后用dlfeval去触发它。这里有个实现细节值得注意。dlgradient要求第一个参数的形状和第二个参数一致或者第一个参数是标量。所以求u对x的偏导时不能直接把X这个2×N的矩阵传进去然后指望它分开算得先把X拆成x和y两个1×N的向量再分别求导。这也是我在下面的代码里先写x X(1, :); y X(2, :);的原因。2.2 用dlnetwork搭建全连接网络网络搭建我用的是dlnetwork这是MATLAB里支持自动微分的深度学习网络容器。网络结构很朴素输入层2个神经元对应(x, y)中间三个隐藏层各50个神经元输出层1个神经元对应u。激活函数必须用tanh不能用ReLU。原因在于PDE损失函数里需要二阶导数而ReLU的二阶导恒等于0。用ReLU做激活函数神经网络的输出对输入的二次导数在定义域内几乎处处为零根本没有能力表征方程里的曲率项。tanh是光滑函数各阶导数都存在而且导数在饱和区会趋近于0但不会完全消失这对梯度传播的稳定性有好处。隐藏层的层数和宽度并非越多越好。我试过单隐藏层和五层隐藏层经验是二维泊松方程这种简单问题上三层隐藏层性价比最高。层数太少网络的拟合能力不足以精确逼近二阶导的复杂组合层数太多小样本情况下容易过拟合训练反而更慢且不稳定。featureInputLayer(2, Normalization, none)这里要特别注意不要开标准化。PINN的输入是物理坐标坐标范围通常很规整加了标准化反而会把坐标映射到不直观的尺度上影响后续对误差的分析。2.3 用adamupdate做优化更新损失函数定义好之后剩下的就是标准的深度学习训练流程。优化器我用的是AdamMATLAB里对应的函数是adamupdate。Adam的优势是自适应学习率对PINN这种多损失项叠加的场景比较友好不需要手动频繁调整学习率。训练策略上有一个经验先跑一万步左右的Adam把损失压到一个比较低的量级然后再根据损失曲线的形态决定是否继续。纯Adam训练在PINN问题上的收敛精度通常能满足基准测试需求但如果要追求更高精度可以考虑切换到L-BFGS这类二阶优化器做精调。MATLAB新版本里lbfgsupdate函数可以直接替代adamupdate实现逻辑几乎不用改。adamupdate维护一组内部状态变量trailingAvg和trailingAvgSq这两个变量要在训练循环外初始化成空数组然后在循环里反复更新。注意学习率不宜设得太大1e-3是一个稳妥的起点我在这篇示例里也是用这个值。3. 一版能直接跑通的完整源码采样、训练、验证3.1 生成训练数据内部点与边界点PINN的数据生成和传统深度学习完全不同不需要任何外部数据集而是在求解区域内随机采样坐标点。内部点用来计算方程残差边界点用来计算边界条件残差。内部点我用rand在[0,1]×[0,1]内均匀采样2000个点。这里用均匀随机采样而不是网格化采样是故意的。网格化采样会造成点的规则分布在某些情况下会让网络学到网格相关的伪模式均匀随机采样则让网络每次都看到略有差异的分布规律泛化更稳。边界点四条边各取100个点一共400个。随机种子设成固定值这一步很重要。PINN训练结果有一定随机性固定随机种子能保证实验结果可复现便于调试和对比。我在代码开头写了rng(42)你们自己调试时也可以固定一个种子。3.2 构建网络与损失函数下面是完整的主程序代码从问题定义到网络搭建我加了详细的注释。%% 一维调试视角下的二维泊松方程PINN求解 % 环境MATLAB R2022a 及以上Deep Learning Toolbox % 问题-u_xx - u_yy 2*pi^2*sin(pi*x)*sin(pi*y) % 区域[0,1]x[0,1]边界 u 0 % 解析解u sin(pi*x)*sin(pi*y) clear; clc; rng(42); %% 1. 问题定义 uExact (x,y) sin(pi*x).*sin(pi*y); fSource (x,y) 2*pi^2*sin(pi*x).*sin(pi*y); %% 2. 生成训练数据 N_in 2000; % 内部点数量 x_in rand(N_in, 1); y_in rand(N_in, 1); X dlarray([x_in; y_in], CB); % 2×NC通道B批量 F dlarray(fSource(x_in, y_in), CB); % 源项值 N_b 100; % 每条边采样点数 xb1 rand(N_b,1); yb1 zeros(N_b,1); % 底边 y0 xb2 rand(N_b,1); yb2 ones(N_b,1); % 顶边 y1 yb3 rand(N_b,1); xb3 zeros(N_b,1); % 左边 x0 yb4 rand(N_b,1); xb4 ones(N_b,1); % 右边 x1 Xb dlarray([xb1; xb2; xb3; xb4], CB); Yb dlarray([yb1; yb2; yb3; yb4], CB); XbAll [Xb; Yb]; % 2×(4*N_b) UB dlarray(zeros(1, 4*N_b), CB); % 边界目标值全零这段代码里有一个容易忽略的点边界点生成时我分别对四条边用rand采样然后拼接成完整边界点集。这样做比在边界上均匀取点更符合PINN的随机采样精神且四条边各自的点数容易控制。如果你想让某条边对解的约束更强可以单独增加那一条边的采样密度。3.3 网络定义与损失函数实现网络定义和损失函数是整段代码的核心。损失函数里涉及二阶导的自动微分计算我在注释里标出了每一步的数学含义。%% 3. 构建神经网络 hiddenSize 50; layers [ featureInputLayer(2, Normalization, none) fullyConnectedLayer(hiddenSize) tanhLayer fullyConnectedLayer(hiddenSize) tanhLayer fullyConnectedLayer(hiddenSize) tanhLayer fullyConnectedLayer(1) ]; dlnet dlnetwork(layers); %% 4. 定义模型损失函数需保存为独立函数文件或嵌套函数 function [loss, pdeLoss, bcLoss] modelLoss(dlnet, X, F, XbAll, UB) % 拆出x和y坐标 x X(1, :); % 1×N 内部点x坐标 y X(2, :); % 1×N 内部点y坐标 % 前向传播得到预测u U forward(dlnet, X); % 一阶导dlgradient要求第一个参数与第二个参数同形状 Ux dlgradient(U, x); Uy dlgradient(U, y); % 二阶导对一阶导再求一次梯度 Uxx dlgradient(Ux, x); Uyy dlgradient(Uy, y); % PDE残差损失-Uxx - Uyy - f pdeResidual -Uxx - Uyy - F; pdeLoss mean(pdeResidual.^2); % 边界条件损失u在边界上应为0 Ub forward(dlnet, XbAll); bcResidual Ub - UB; bcLoss mean(bcResidual.^2); % 总损失 loss pdeLoss bcLoss; end这里为什么可以对一阶导Ux再次调用dlgradient这是MATLAB自动微分能力的体现。Ux本身是从计算图里求出来的dlarray它保留了前向传播的所有运算轨迹。再一次调用dlgradient(Ux, x)时MATLAB会在原来的计算图上继续做反向传播从而得到二阶导数。理解这一点就理解了PINN在MATLAB里的全部技术核心。关于损失权重这里PDE残差和边界残差用的是默认的1:1权重。对这个问题零Dirichlet边界且边界值精确等于解析解在边界上的限制1:1能收敛得很好。如果遇到边界损失显著大于PDE损失的情况可以考虑给边界损失一个大于1的权重系数我后面会详细讲。3.4 训练主循环与训练过程监控训练循环用Adam优化器迭代20000步每1000步打印一次损失值。%% 5. 训练主循环在主脚本中执行 numIter 20000; learnRate 1e-3; trailingAvg []; trailingAvgSq []; lossHistory zeros(numIter, 1); pdeLossHistory zeros(numIter, 1); bcLossHistory zeros(numIter, 1); for iter 1:numIter % dlfeval触发带自动微分的损失计算和梯度反传 [loss, grads, pdeLoss, bcLoss] dlfeval(modelLoss, dlnet, X, F, XbAll, UB); % Adam更新网络参数 [dlnet, trailingAvg, trailingAvgSq] adamupdate(dlnet, grads, ... trailingAvg, trailingAvgSq, iter, learnRate); lossHistory(iter) extractdata(loss); pdeLossHistory(iter) extractdata(pdeLoss); bcLossHistory(iter) extractdata(bcLoss); if mod(iter, 1000) 0 fprintf(Iter %6d | Loss: %.3e | PDE: %.3e | BC: %.3e\n, ... iter, lossHistory(iter), pdeLossHistory(iter), bcLossHistory(iter)); end endextractdata的作用是把dlarray里的数值取出来变成普通的double这样才能在循环里做数组存储和画图。打印的信息里我同时输出了PDE损失和BC损失这两个值的相对大小很值得观察它能告诉你训练是收敛到了平衡状态还是某一项损失被过度压制了。我实际跑出来的损失曲线是前几百步损失从初始的O(1)量级迅速下降到1e-3左右之后下降速度放缓到两万步时总损失在1e-5量级。PDE损失和BC损失大致同步下降没有出现某一项长期压制另一项的情况。4. 在验证集上看效果误差分布与收敛性4.1 预测解、解析解、绝对误差三维图对比训练完成后需要验证网络学到的函数到底准不准。验证点我用网格来生成为了方便画三维曲面图取0.02为步长在这个密集网格上对比预测值和解析解。%% 6. 验证与可视化 [Xg, Yg] meshgrid(0:0.02:1, 0:0.02:1); Xgv Xg(:); Ygv Yg(:); % 网格点上的预测解 U_pred predict(dlnet, dlarray([Xgv; Ygv], CB)); U_pred reshape(extractdata(U_pred), size(Xg)); % 解析解 U_true uExact(Xg, Yg); % 绝对误差 err abs(U_pred - U_true); % L2相对误差 L2rel norm(U_pred(:) - U_true(:)) / norm(U_true(:)); fprintf(L2 相对误差: %.4e\n, L2rel); % 最大误差 maxErr max(err(:)); fprintf(最大绝对误差: %.4e\n, maxErr); %% 7. 画图 figure(Position, [100 100 1200 350]); subplot(1, 3, 1); surf(Xg, Yg, U_pred, EdgeColor, none); title(PINN 预测解); xlabel(x); ylabel(y); zlabel(u); subplot(1, 3, 2); surf(Xg, Yg, U_true, EdgeColor, none); title(解析解); xlabel(x); ylabel(y); zlabel(u); subplot(1, 3, 3); surf(Xg, Yg, err, EdgeColor, none); title(绝对误差); xlabel(x); ylabel(y); zlabel(error);我用这套代码跑出来的典型结果是L2相对误差在1e-3到1e-4之间最大绝对误差出现在边界附近。这个精度虽然没有有限差分法在精细网格上的精度高但已经足以验证框架和代码的正确性。误差云图呈现出的规律一般是边界处误差稍大中心区域误差更小误差大的区域往往也是解变化剧烈的区域。注意这里用的是predict而不是forward。predict只做前向传播不构建计算图推理速度快且不占用自动微分的内存开销。验证阶段不需要求导所以用predict是标准做法。4.2 损失曲线解读什么时候算训练好了训练过程中保存的损失曲线是判断收敛状态的重要依据。我把损失曲线画出来的经验是这样的理想情况下PDE损失和BC损失同时下降最终都稳定在一个较小的量级。如果总损失下降但PDE损失长期高于BC损失很多说明网络优先去迎合边界约束但内部的方程满足度不够。这时候需要增加内部采样点的数量或者给PDE损失一个更大的权重。如果BC损失长期高于PDE损失说明边界条件没有被充分学习一般是因为边界点太少或边界损失权重太低。“训练好了”的判断标准要看你最终的应用场景。如果只是验证算法可行性损失降到1e-4就完全够了。如果要做高精度计算那么损失量级要更低并且最好用独立验证集的误差来确认而不能只看训练损失。5. 从基准问题走向实际应用调试清单与扩展方向5.1 五条最值得记住的实操经验我把跑通这个完整流程之后踩过的坑和积累的经验整理成了一张表先看现象再给处理方案。问题现象根本原因处理方式训练NaN学习率过大或参数初始化范围过大学习率降到5e-4以下或者改用initialize重新初始化网络损失降到1e-3就卡住网络容量不足或采样点太少增加隐藏层宽度到80或增加到4000个内部点PDE损失和BC损失数量级差很大两类损失权重失衡在总损失中给较小的一项乘以大于1的权重系数边界附近误差特别大边界条件约束不够强增加边界采样点数量或提高边界损失权重多次训练结果差异大PINN随机性大未固定种子或训练不充分固定随机种子增加训练迭代次数多跑几次取最优这里有一个经验值得展开固定随机种子。PINN的求解结果对初始化比较敏感同样的代码换一个随机种子最终L2误差可能相差一个数量级。所以我的习惯是正式实验固定rng(42)但如果要判断某个超参数改动的真实效果会跑3到5个不同的种子取平均避免被单次随机性误导。5.2 从二维泊松方程扩展出去的几条路这个基准问题跑通之后扩展方向其实非常多我梳理了四条我认为最容易上手且工程价值最高的路径。第一改变源项和边界条件类型。从Dirichlet边界换成Neumann边界损失函数里需要加一项边界法向导数残差从零边界换成本质边界只需要把UB从零向量改成对应的边界函数值。这是PINN最舒服的改动方式代码层面的变化很小。第二不规则区域。这是PINN相对传统方法最有优势的地方。区域改成L形或者圆形不需要重新剖网格只需要把内部采样点的采样范围限制在区域内。用inpolygon之类的函数做点在区域内的判断就能生成任意形状区域的训练点。第三带时间项的抛物型方程。在输入维度上增加时间t网络输入从(x,y)变成(x,y,t)损失函数在PDE残差里多一项时间导数。很多做瞬态热传导和反应扩散方程的人都用这个思路。第四参数化PDE。把方程中的某个系数比如扩散系数也作为网络输入的一部分训练好的网络可以在给定系数后直接输出对应解相当于一次训练、多参数求解。我自己的实际体会是从基准问题扩展的第一步最推荐做的是换一个非零Dirichlet边界比如边界上u(x,0)x这种简单的线性分布。这个改动会让你深刻理解边界残差在损失函数里的作用方式比直接上复杂问题要顺利得多。5.3 我对PINN选型的一点个人判断围绕PINN本身我需要说句公道话。它并不是要取代传统数值方法而是一种互补的工具。传统有限元、有限差分在规则区域、高精度需求下仍然是无可争议的首选。PINN的价值在于处理传统方法麻烦的场景不规则区域、反问题根据部分观测反推方程系数、参数化工况族以及和实验数据直接结合的物理-数据混合建模。对于刚接触PINN的读者我的建议是先把这篇文章的基准问题完整跑通加入自己的注释手动改几个超参数观察效果然后从这个稳定起点出发去扩展。不要一上来就挑战复杂工程问题先把四个核心环节——采样、网络、损失、优化——的直觉建立起来。后续在实际应用里需要在这个基础上叠加什么技巧注意力权重、自适应采样、多尺度网络也就有稳固的根基了。这套MATLAB源码我在自己的测试机上跑了多个版本R2022a、R2023b都能稳定运行。如果你用的是老版本注意adamupdate是R2019b之后才有的函数太老的版本需要换成手动实现Adam更新公式。最后再分享一个细节训练完成后把网络参数保存下来后续预测时直接load进来用不用重新训练这才是真正把模型用起来的姿势。本文还有配套的精品资源点击获取
返回列表