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、dlfeval

MATLAB做PINN，核心工具是Deep Learning Toolbox里的自动微分能力。三个关键函数要搞明白：dlarray用来封装带计算图的数据，dlgradient用来求偏导，dlfeval用来在计算图环境下执行函数。

PINN的损失函数里需要算u对x和y的二阶偏导。神经网络的前向传播本身是一个复杂的复合函数，手动推导二阶导完全不现实。自动微分会在前向传播时记录每一步运算关系，然后通过链式法则自动求出各阶导数。在MATLAB里，这个过程被抽象成了dlgradient的两次调用：第一次对前向传播结果求一阶导，第二次对一阶导再次求导得到二阶导。

一个容易踩的坑是：dlgradient不能直接在普通函数里调用，必须在dlfeval的上下文里执行。我第一次写的时候，在普通脚本里直接调用dlgradient，MATLAB直接报错，提示只能在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做优化更新

损失函数定义好之后，剩下的就是标准的深度学习训练流程。优化器我用的是Adam，MATLAB里对应的函数是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×N，C=通道，B=批量 F = dlarray(fSource(x_in, y_in)', 'CB'); % 源项值 N_b = 100; % 每条边采样点数 xb1 = rand(N_b,1); yb1 = zeros(N_b,1); % 底边 y=0 xb2 = rand(N_b,1); yb2 = ones(N_b,1); % 顶边 y=1 yb3 = rand(N_b,1); xb3 = zeros(N_b,1); % 左边 x=0 yb4 = rand(N_b,1); xb4 = ones(N_b,1); % 右边 x=1 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 end

extractdata的作用是把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 五条最值得记住的实操经验

我把跑通这个完整流程之后踩过的坑和积累的经验整理成了一张表，先看现象，再给处理方案。

这里有一个经验值得展开：固定随机种子。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进来用，不用重新训练，这才是真正把模型用起来的姿势。

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