这个项目表面上只是把几个词拼在一起——傅里叶特征、PINN、一维Burgers方程、MATLAB代码——实际上它是一条完整的深度学习求解PDE的路径。我最初就是被这六个词带进坑的写出来的网络怎么调也收敛不到参考解等真正搞明白傅里叶特征和损失函数配比之后才发现PINN的工程实现比算法本身更容易踩坑。这篇文章整理的是我用MATLAB从零搭完这整套代码的完整过程包括损失函数怎么拆、自动微分怎么用、频率怎么选、最后怎么判断结果靠不靠谱。适合已经有深度学习基础、想把PINN落到实际代码里的读者也适合被MATLAB自动微分和各种报错折磨的人。1. 问题背景与整体设计思路1.1 为什么拿一维Burgers方程当测试床一维Burgers方程长这样u_t u * u_x ν * u_xx其中 ν 是黏性系数u 是速度场x 是空间坐标t 是时间。它同时包含非线性对流项u*u_x和二阶扩散项ν*u_xx既能模拟出激波一样的陡峭梯度又不像Navier-Stokes方程那样涉及速度和压力的耦合是检验数值算法最合适的“训练场”。我做这版代码时使用了经典参数ν 0.01/πx 取 [-1, 1]t 取 [0, 1]。初始条件取u(0,x) -sin(π*x)边界条件为u(t,-1) u(t,1) 0。这个设定之所以在PINN论文里反复出现是因为它在 t 约等于 0.35 之后会在 x 0 附近形成一个很陡的间断层网络要想把这个间断层拟合出来就不得不同时处理低频平滑区域和高频陡峭区域这对标准神经网络来说相当不友好。如果只用有限差分或者有限元去解这个方程网格加密是主要手段但PINN的思路完全不同它把解 u(x,t) 直接参数化为一个神经网络输入是坐标 (x,t)输出是速度值 u然后用自动微分计算方程里的偏导数通过最小化方程残差来训练网络。这样做的好处是不需要网格划分也不需要构造离散格式适合更复杂的几何区域代价是对训练技巧要求很高。1.2 为什么标准全连接网络在PDE任务上不够用很多人第一次写PINN都会犯一个同样的错误以为只要网络足够宽、训练足够久就能逼近任意函数。对简单的Poisson方程可能确实如此但一旦遇到Burgers方程这种带激波的问题普通全连接网络会表现得非常“软”激波附近的梯度被抹得又宽又平怎么加大网络容量都改善有限。这背后的核心问题是“频谱偏差”。全连接网络在实际训练中会优先学到低频分量高频分量要耗费大量迭代次数才会慢慢出现。而Burgers方程的激波区域恰恰是高频信息最集中的地方。模型在初期只拟合大尺度轮廓后期却已经收敛得过早激波位置的高频细节被卡死了。解决办法就是给输入坐标加“傅里叶特征”。简单说不直接把 (x,t) 喂给网络而是先把它们映射到一组不同频率的三角函数上cos(2πf x), sin(2πf x), cos(2πf t), sin(2πf t), ...这样网络在输入层就已经拿到不同频段的信息梯度更新的过程中低频和高频分量可以被同时驱动起来。视觉上可以把它理解成给坐标加了一组“频率坐标指纹”网络不再需要自己硬生生地从零构造频率分解。这个思路最开始在NeRF里被广泛使用后来被引入PINN后对激波类问题的改善非常明显。1.3 为什么选MATLAB而不是Python现在PINN生态里Python的工具链是最全的PyTorch、TensorFlow、DeepXDE各有优势但我这版项目最终选了MATLAB有几个实际原因。首先是MATLAB深度学习工具箱原生支持dlarray和dlgradient自动微分不需要额外装一套神经网络框架代码结构非常干净适合快速验证。其次是MATLAB的调试体验好我在写损失函数和梯度传播时可以在任意一行停下来查看张量形状这对排查二阶导数的维度问题帮助很大。另外如果后续要把这个求解器接到控制系统或信号处理任务上MATLAB/Simulink的天然生态会让集成容易很多。我见过不少做电池建模、传感器校准的工程师他们手里就是一套Simulink模型这时候让他为了跑个PINN再单独搭一套Python环境维护负担太重。所以选择MATLAB不是因为它比Python强而是它更适合具体场景里“快速跑通、接上流程”的需求。2. 关键原理损失函数、自动微分与傅里叶特征2.1 损失函数的具体构成PINN的核心不是网络结构有多复杂而是损失函数怎么搭。Burgers方程没有现成的监督标签我们能利用的是物理约束内部点满足方程初始点满足初始条件边界点满足边界条件。所以损失函数由三部分组成L L_residual L_initial L_boundaryL_residual 是方程残差损失对内部采样点 (x_r, t_r) 计算r u_t u * u_x - ν * u_xx L_residual mean(r^2)L_initial 是初始条件损失L_initial mean( ( u(0, x_0) sin(π*x_0) )^2 )L_boundary 是边界条件损失L_boundary mean( u(t, -1)^2 u(t, 1)^2 )默认情况下三者权重都取1但我在实验里发现如果初始条件和边界条件的采样点太少损失会被残差项压住网络会“无视”边界。最简单有效的做法是先保证初始和边界采样点数达到几百个让三项损失在数量级上能互相被看到。内部残差点的采样方式也很讲究。我一开始用均匀网格采样发现激波区域拟合效果一般改成随机采样后训练过程中每次采到的点分布都不同相当于不断给网络提供新的局部约束收敛稳定性明显更好。后面如果想再进一步可以用残差自适应采样把采样点往损失大的区域集中。2.2 MATLAB中的自动微分与二阶导MATLAB的自动微分是围绕dlarray设计的。你需要把输入数据包装成dlarray然后所有对网络参数和输入的计算都会生成一个计算图通过dlgradient得到导数并且必须在dlfeval里执行。对一维Burgers方程来说我们不仅需要一阶导数u_t和u_x还需要二阶导数u_xx这会带来一个小坑。dlgradient只计算某个输出对输入的一阶梯度。想要二阶导数必须嵌套调用也就是先对一阶导数再求一次梯度。我在代码里是这样处理的U forwardNet(X, params, freqs); % 一阶导数 grads dlgradient(sum(U, all), X); Ux grads(1, :); Ut grads(2, :); % 二阶导数 Ux_temp dlgradient(sum(Ux, all), X(1, :)); Uxx Ux_temp(1, :);这里最关键的一点是每次dlgradient都要求目标值是标量所以用sum(..., all)把所有样本加成一个标量。一开始我直接对向量做dlgradient(U, X)MATLAB直接报错后来查文档才发现自动微分接口的约束。这个细节很容易被网上的教程忽略但在MATLAB里是绕不过去的。另外注意嵌套求导会让计算图变深内存占用比单次求导要高。内部残差点数量设到一万点以后GPU显存不够的话可以适当降到五千点训练效果差别不会太大但内存压力能缓解不少。2.3 傅里叶特征的频率选择傅里叶特征不是越密越好频率选多少是个需要权衡的问题。我常用的特征构造函数长这样function psi fourierFeature(X, freqs) % X: 2xN dlarray第一行是x第二行是t psi []; for k 1:numel(freqs) f freqs(k); psi [psi; cos(2*pi*f*X); sin(2*pi*f*X)]; end psi [X; psi]; endfreqs我习惯用2.^(0:N-1)生成也就是频率呈倍数递增。比如N6时得到[1, 2, 4, 8, 16, 32]。这样可以覆盖从低频到高频的多个尺度等于给网络提供了一组位置的频率描述。实际测试下来频率太少比如只有1个频率对激波没有明显改善频率太多比如到128以上会导致网络在训练初期疯狂震荡损失函数像心电图一样上下跳收敛非常困难。我最后常用的区间是4到8个频率最大频率控制在32左右。这个范围对一维Burgers方程已经足够再往上加就是纯增加计算量。3. 完整MATLAB实现与代码解析3.1 工程文件组织整个项目我拆成了四个文件方便在后续换方程时复用文件作用main_burgers.m主脚本负责参数设置、数据采样、训练network.m网络前向计算含傅里叶特征映射model_loss.m损失函数返回总损失和梯度plot_results.m后处理与可视化主脚本负责串起所有流程network.m只管前向计算model_loss.m被dlfeval调用并返回梯度这样的结构比把所有代码堆在一个脚本里调试起来好很多。3.2 网络初始化与前向传播网络结构很常规输入层接傅里叶特征中间四层全连接隐藏层使用tanh激活输出层是线性层。隐藏层神经元数量我设为20这个规模应对一维问题已经足够。初始化权重时不要全部用同一个分布最好使用Xavier/Glorot初始化让每一层的输入输出方差保持一致。我写成了这样% network.m function U network(X, params, freqs) H fourierFeature(X, freqs); % 输入层加傅里叶特征 for k 1:numel(params.W)-1 H tanh(params.W{k} * H params.b{k}); end U params.W{end} * H params.b{end}; end function params initNetwork(layerDim) params.W {}; params.b {}; for k 1:numel(layerDim)-1 fanIn layerDim(k); fanOut layerDim(k1); sigma sqrt(2 / (fanIn fanOut)); % Xavier params.W{k} dlarray(randn(fanOut, fanIn, single) * sigma); params.b{k} dlarray(zeros(fanOut, 1, single)); end end激活函数我固定用tanh不用ReLU。原因是ReLU在负半轴的梯度恒为0迭代过程一旦某个神经元落入负区间它会长期“休眠”对PDE这种需要平滑梯度的拟合任务非常不利。tanh处处可导且关于原点对称在残差损失传播时更稳定。如果想追求更高精度可以把tanh换成sin或swish变体但默认tanh是性价比最高的选择。3.3 损失函数与训练循环损失函数部分的代码量比较大核心就是把三个损失项分别算出来再求和。我做成了独立函数function [loss, gradients] model_loss(params, X_res, T_res, X_ic, U_ic, X_bc, U_bc, nu, freqs) % 内部残差点的loss U network([X_res; T_res], params, freqs); dU dlgradient(sum(U, all), [X_res; T_res]); U_x dU(1, :); U_t dU(2, :); dU_x dlgradient(sum(U_x, all), X_res); U_xx dU_x(1, :); residual U_t U .* U_x - nu * U_xx; loss_res mean(residual.^2, all); % 初始条件loss U_ic_pred network([X_ic; zeros(size(X_ic))], params, freqs); loss_ic mean((U_ic_pred - U_ic).^2, all); % 边界条件loss U_bc_pred network([X_bc; T_bc], params, freqs); loss_bc mean((U_bc_pred - U_bc).^2, all); loss loss_res loss_ic loss_bc; gradients dlgradient(loss, params); end这里有个细节网络输入是把x和t拼接成一个[2, N]的张量第一行是空间坐标第二行是时间坐标。这样做是因为网络需要同时接收两个变量并且后续对它们分别求偏导时更方便。如果你习惯把坐标拼成一维向量也行但维度处理会更绕不建议写成那样。训练循环用Adam优化器MATLAB里直接调用adamupdate就行for iter 1:numEpochs [loss, grads] dlfeval(model_loss, params, X_res, T_res, ... X_ic, U_ic, X_bc, U_bc, nu, freqs); learnRate initialLearnRate * (0.5 ^ floor(iter / 2000)); [params, avg, avgSq] adamupdate(params, grads, avg, avgSq, iter, learnRate); end学习率我从1e-3起步每2000轮衰减一半总共跑6000到8000轮就基本收敛。这个衰减策略比固定学习率稳定因为训练初期需要一个相对大的步长快速进入有利区域后期则需要更小的步长去精修激波位置。如果损失曲线在中途出现明显的周期性波动说明学习率偏大要适当调低初值或者把衰减周期调短。3.4 结果可视化与误差统计训练完以后光看最终损失值不够必须把预测解和参考解放到一起对比。我通常在三个时间截面t 0.25, 0.5, 0.75上把预测结果画出来并计算L2相对误差function plot_results(X_grid, T_grid, U_pred, U_ref, errL2) figure; for k 1:3 tVal [0.25, 0.5, 0.75](k); idx abs(T_grid - tVal) 1e-6; plot(X_grid(idx), U_pred(idx), r-, X_grid(idx), U_ref(idx), b--); end legend(PINN, Reference); xlabel(x); ylabel(u); title(sprintf(L2 error: %.3e, errL2)); end参考解怎么来最稳的方法是先用一个高分辨率的有限差分法算一遍作为基准或者直接用MATLAB自带偏微分方程工具箱求一个近似解。我习惯直接对网格点用伪谱法或高精度差分生成参考值然后在验证点上做插值。这样算出来的L2误差能客观反映PINN的精度一般在0.01以下就算表现不错了。4. 训练配置、参数调节与结果经验4.1 一组可以直接开跑的默认参数参数取值说明ν0.01/πBurgers方程黏性系数空间范围[-1, 1]一维空间域时间范围[0, 1]时间域傅里叶频率数62.^(0:5)隐藏层4层×20神经元tanh激活内部残差采样点8000随机均匀分布初始条件采样点200均匀采样边界条件采样点200每条边界100个优化器Adamadamupdate初始学习率1e-3每2000轮衰减0.5训练轮数8000视损失下降情况调整这套参数我实测下来能在几分钟内跑完CPU训练用GPU会更快。关键在于残差采样点不是越多越好训练初期用8000个点已经能让损失下降到1e-4量级再增加到两万点只会拖慢每轮迭代精度提升却非常有限。4.2 如何判断训练是否正常训练过程中我会同时盯三个东西总损失曲线、残差损失曲线、验证点L2误差。一个好的趋势是总损失从初始的0.1量级平滑下降到4000轮左右降到1e-3附近之后进入慢速精修阶段L2误差随损失一起下降。如果只盯总损失很容易被骗。我遇到过损失已经降到1e-5但预测曲线在激波处明显偏离参考解的情况。原因在于激波区域宽度很窄占整个求解域的样本比例很小均方误差被平滑区域主导激波附近的局部误差在总损失里体现不出来。所以必须隔一段时间就在验证网格上算一次真实误差。另外初始条件损失如果下降得比残差损失慢网络会更倾向于“违背”初始条件最终结果在早期时刻会整体偏移。这时可以把初始和边界损失权重从1调到5或10强制网络优先满足这两个条件。4.3 频率和网络宽度的影响配置现象结论无傅里叶特征激波被抹平L2误差在0.05以上高频信息学不到位频率数3激波位置基本正确但峰值偏低频率覆盖不够密频率数6激波和峰值都能对上误差最优推荐配置频率数10前期训练震荡收敛变慢频率过高适得其反隐藏层20个神经元够用一维问题不需要过宽隐藏层80个神经元训练变慢但精度没有明显提升过拟合风险上升网络宽度不是越宽越好。PINN把方程约束放进损失函数后模型容量过大会让网络更容易把残差损失“死记硬背”到某个局部最小值泛化反而变差。我最后固定用20到40个神经元在保证精度的同时减少调试成本。5. 常见问题与排查实录5.1 MATLAB中文注释乱码这个问题最容易让人抓狂从网上下载或用旧版本保存的.m文件中文注释在新版MATLAB里变成一串乱码代码本身没问题但看着就脑壳疼。原因基本是文件编码不匹配。老版本MATLAB在中文系统里默认用GBK/ANSI保存.m文件而新版MATLAB默认用UTF-8打开编码对不上自然乱码。解决方法有两个。第一个是打开文件时手动指定编码open(myfile.m, encoding, UTF-8)如果文件本身是GBK编码就指定GBK。第二个是改编辑器默认编码在顶部菜单找到“预设”或“Preferences”进入“MATLAB 编辑器/调试器 语言”把文件编码改成“UTF-8”然后重新打开文件。对已经在磁盘上乱码的文件我一般先用系统自带记事本打开另存为UTF-8编码再回到MATLAB打开基本都能恢复。5.2 训练出NaN、梯度爆炸训练一两千轮后损失突然变成NaN最常见的原因是学习率偏大或网络初始化范围太大。我在初始化时把权重标准差控制在一定范围内后NaN出现频率大幅下降。另一个隐蔽原因是残差损失里出现了无穷大梯度尤其在激波刚开始形成、梯度变化非常剧烈的阶段。排查时先做三件事把学习率降到1e-4重新初始化网络并把内部采样点重新随机生成。如果NaN还在就要检查傅里叶特征的最大频率是否过高。频率超过64后正弦和余弦的梯度变化极快网络很难稳定追踪这时把最大频率降到32以下基本能解决问题。5.3 损失下不去或预测偏差如果损失始终停在0.1量级不下降先怀疑是损失权重失衡。内部残差项的计算量远大于初始和边界项如果后者采样点太少它们的梯度很容易在反向传播中被淹没。解决方案是把初始和边界损失权重调大或者增加采样点数量让三种约束在梯度层面达到平衡。如果损失已经很低但预测和参考解对不上那就不是优化问题而是“验证方式”问题。可能你在训练时用了固定网格采点而激波恰好落在两排采样点之间的空隙里。改成随机采样并每隔一段时间重新生成采样点能明显缓解这种局部遗漏。5.4 自动微分相关报错自查MATLAB的自动微分报错几乎都在我刚接触的时候出现过最常见的几个场景是dlgradient被调用在dlfeval之外报错说无法解析梯度。检查一下代码结构确保所有自动微分都包在dlfeval里面。目标值不是标量报错说梯度形状不匹配。解决办法是sum(U, all)把结果压成标量再求导。输入不是dlarray报错说变量未启用自动微分。检查输入是否已经用dlarray()包装过。输入维度对不上比如params.W{k}和H的矩阵乘法维度不一致。建议在网络前向函数入口加一行size(H)调试输出被报错折磨时这招最有效。6. 扩展方向与个人实操体会6.1 从一维到更高维这套代码从一维Burgers方程扩展到二维Burgers方程或者稳态Navier-Stokes方程核心改动很小。主要是把输入从(x,t)变成(x,y,t)并增加u、v两个速度分量对应的残差方程。网络输出层维度从1变成2或3损失函数里多拼几个残差项其它训练机制完全不用动。6.2 训练策略可以继续升级如果追求更高精度可以在现有基础上加入残差自适应采样每训练一定轮数后重新采样把残差大的区域作为新的采样点。另一个思路是使用两阶段训练先用Adam跑一个粗略解再切到L-BFGS类优化器做精细优化。后者在MATLAB里需要自己实现或借助外部工具稍微麻烦但对精度提升很直接。6.3 个人体会PINN这项目最容易被低估的地方不是算法本身而是工程细节。频率怎么编码、损失权重怎么配、采样点怎么分布、自动微分怎么保证标量输出每一个细节都能让结果从“跑得起来”变成“真的准确”。我把这套代码在MATLAB上跑通之后最大的收获是明白了物理约束不该只被当作一个正则项它其实是在告诉网络“哪些函数形式不可能出现”这种约束比任何数据增强都来得强。GitHub上类似的Python实现很多但MATLAB版本相对少所以我一直在维护这套代码后续我会继续补充二维问题和不同参数下的实验对比如果你也在做类似的工作可以对照本文的代码框架做修改欢迎在评论区交流你的训练结果。