一维Burgers方程这个题目相信做流场仿真或者偏微分方程数值解的朋友都不陌生。它式子很简洁只有非线性对流项和粘性项但解里能长出陡峭的激波是测试数值算法最经典的“试金石”。这两年物理信息神经网络PINN几乎成了用AI解PDE的代名词顶刊、开源库、公众号里到处都是它的身影可公开实现大多是Python版很多熟悉MATLAB的工程师和研究生想上手还得先临时转语言、配环境。这篇博文就把我调好的一份基于傅里叶特征Fourier Feature的PINN求解一维Burgers方程的MATLAB代码完整拆开讲一遍想直接跑通的人可以按步骤抄作业想弄明白原理的人也能从损失函数设计、Fourier特征映射、自动微分求残差这些关键点看清每一步背后到底在干什么。1. 项目概述与整体设计思路1.1 这个项目到底做了什么项目的目标可以浓缩成一句话用神经网络去逼近一个二维函数 u(x,t)让这个函数在定义域内满足Burgers方程同时在初始时刻满足 -sin(πx)在左右边界满足零值约束。与常规数值方法不同的是这里没有网格、没有差分模板、没有有限元刚度矩阵所有约束都以损失项的形式加入训练过程网络每更新一次参数就是在“解一次方程”。这种思路的直观含义是神经网络不直接学习数据标签而是学习“物理规律本身”。内部配点处方程残差越小说明网络的输出越接近真实解初边值点的误差越小说明解与定解条件贴合越好。最终训练完成后网络就成为一个可微的函数逼近器输入任意 (x,t)输出对应的速度 u。1.2 为什么选一维Burgers方程做基准选Burgers方程不是因为简单而是因为它“简单却难点密集”。方程左边是经典的非线性对流项右边是粘性扩散项两者竞争会让解在某个时刻之后出现接近间断的陡峭激波。在 ν0.01/π 这种低粘性设置下激波会非常锋利普通数值方法处理起来都相当考验更不用说神经网络。用这个方程还有两个现实原因。第一它有相对成熟的参考解获取途径Cole-Hopf变换可以得到解析解MATLAB的 pdepe 也可以算高精度参考解验证误差很方便。第二Raissi在2019年提出PINN时的原始benchmark就是它各种改进方法都会拿这个方程做横向对比所以我这套加Fourier Feature的实现在论文里也能找到对应参照。1.3 为什么必须加Fourier Feature这是整个项目最关键的设计决定不是锦上添花。普通全连接网络在逼近函数时存在明显的低频偏置也就是从优化理论的视角看网络在训练初期只会优先拟合低频成分高频分量更新非常缓慢。Burgers方程的激波恰恰是高频分量集中的区域如果直接用朴素MLP训练结果往往在激波处一片光滑怎么加网络宽度和深度都改善有限。Fourier Feature的做法非常直接输入坐标 (x,t) 先经过一层固定权重的随机余弦/正弦映射将低频坐标变换成多频率组合的特征向量再送入网络。打个比方这相当于给相机同时装上了广角镜头和长焦镜头——网络既能处理平滑区域的低频变化也有能力刻画激波附近的高频细节。实测下来加上这层之后激波附近的误差能明显下降收敛速度也更快。1.4 为什么用MATLAB而不是Python当前PINN生态确实是Python占主导但MATLAB有一个被低估的优势自动微分和自定义训练循环已经被封装得很顺手了。dlnetwork定义网络、dlgradient求各阶偏导、adamupdate做参数更新三件套配合得当。对于熟悉Simulink、FDTD、有限元工具箱的工程研究人员保留在自己的工作环境里完成神经网络求解不用为了一个算例去安装Python环境、管理依赖包。当然代价也很明显网上针对“MATLAB PINN”的完整教程少得可怜很多细节需要自己踩坑。这也是我写这篇博文的核心动机把可运行的代码、需要避开的坑、以及调参的逻辑一次性讲明白。2. 理论拆解方程、损失函数与傅里叶特征原理2.1 控制方程与定解条件一维带粘性Burgers方程的标准形式为$$\frac{\partial u}{\partial t} u \frac{\partial u}{\partial x} \nu \frac{\partial^2 u}{\partial x^2}$$其中 u(x,t) 是速度场ν 是粘性系数。本项目采用经典算例设置x ∈ [-1,1]t ∈ [0,1]初值 u(x,0) -sin(πx)边界条件 u(-1,t) u(1,t) 0粘性系数 ν 0.01/π。从物理图像上理解非线性对流项 u·u_x 会让波形在传播过程中不断“自陡化”而粘性项尝试把梯度拉平两个过程竞争的结果就是在 t≈0.5 后出现一个非常陡峭的激波且位置大约在 x0 附近。这个解的空间频率跨度极大恰好用来检验神经网络在低频平滑区和高频陡峭区之间的平衡能力。2.2 PINN损失函数的三件套PINN的训练目标是让输出函数在定义域内满足方程在边界初值处满足定解条件。所以损失函数由三部分组成。一是内部配点上的PDE残差把网络输出代入原方程残差为 r u_t u·u_x - ν·u_xx损失取所有配点上残差平方的均值。二是初始条件损失在 t0 的初始点上比较网络输出与 -sin(πx) 的差异。三是边界条件损失在 x-1 和 x1 的边界点上要求网络输出接近0。总损失简单加和即可。三个部分地位可以理解为“三类教师”PDE残差强迫网络遵守物理规律初值项告诉网络起点状态边界项钉住区域边缘。缺任何一项问题都会退化成不适定问题。训练时最大的优势是方程偏导全部交给自动微分用两次 dlgradient 就能拿到一阶和二阶空间导数完全不需要手写差分模板。2.3 Fourier Feature的数学原理Fourier Feature的思想来自NeRF等位置编码工作后来被证明能有效缓解神经网络在低维坐标输入上的频谱偏置问题。映射形式为$$\gamma(v) [\cos(2\pi Bv), \sin(2\pi Bv)]$$其中 v(x,t) 是二维坐标B 是一个 m×2 的随机矩阵每个元素从高斯分布 N(0,σ²) 中采样m 是特征维度的一半。这样每个输入坐标都被展开成一组覆盖不同频率的余弦和正弦基网络在原始坐标下学不到的高频分量在特征空间里从一开始就是显式存在的。关键参数是 σ它决定了频率带宽。σ 太小特征基本只覆盖低频等于白加σ 太大高频噪声过多训练容易震荡甚至发散。按照经验σ 在 5 到 10 这个量级对Burgers方程的激波问题比较合适。至于B矩阵为什么固定而不是训练中更新因为固定随机特征已经可以给网络提供足够的频率工具省掉可学习参数的额外复杂度。2.4 参数选择与网络结构网络结构我采用4层全连接隐藏层各50个神经元激活函数选tanh而不是ReLU。原因很简单输出是要满足二阶导数的ReLU的导数不平滑二阶信息会直接断掉tanh处处光滑能稳定提供 u_xx。Fourier特征维度 m_fourier100所以输入层神经元数是2×100200。配点方面内部配点10000个、初始条件点500个、单侧边界点200个。这个数量在入门配置下足以稳定收敛也不会让训练时间长得让人失去耐心。优化器用Adam初始学习率1e-3迭代8000次。如果用更激进的配置比如加大σ、加深网络可以适当降低学习率防止震荡。参数汇总如下参数取值选择理由ν0.01/π经典基准配置激波足够陡m_fourier100特征维度适中输入200维σ5频率带宽可覆盖激波高频隐藏层3×50深度和宽度平衡激活函数tanh二阶导光滑稳定配点数10000/500/200内部PDE占主导初边值点足够优化器Adam入门稳定无需额外配置学习率1e-3与网络规模和损失量级匹配迭代8000普通桌面CPU可接受3. MATLAB实操从网络结构到训练循环废话少说代码为王。下面直接贴核心代码并逐段说明。整套代码可以按五个块理解参数设置、采样、网络定义、训练、可视化。运行前请先确认已经安装了Deep Learning Toolboxdlnetwork、dlgradient、adamupdate这些函数都依赖它。3.1 参数设置、Fourier特征矩阵与网络构建代码里我把Fourier的2π直接合并进了B矩阵这样后面映射函数就可以少写一个系数。B矩阵用固定随机种子生成保证每次跑结果可复现。网络用featureInputLayer接三层全连接输入特征维度正好是2×m_fourier。%% 1. 参数设置与网络定义 clear; clc; rng(42); % 方程与算例参数 nu 0.01 / pi; % 粘性系数 m_fourier 100; % Fourier 特征维度的一半 sigma 5; % B 矩阵标准差控制频率带宽 Nf 10000; % 内部配点数 Nic 500; % 初始条件采样点数 Nbc 200; % 单侧边界采样点数 nIter 8000; % Adam 迭代次数 lr 1e-3; % 学习率 % Fourier 特征矩阵 B尺寸 m_fourier x 2 B randn(m_fourier, 2) * (2 * pi * sigma); % 构建全连接网络 layers [ featureInputLayer(2 * m_fourier, Normalization, none, Name, in) fullyConnectedLayer(50, Name, fc1) tanhLayer(Name, tanh1) fullyConnectedLayer(50, Name, fc2) tanhLayer(Name, tanh2) fullyConnectedLayer(50, Name, fc3) tanhLayer(Name, tanh3) fullyConnectedLayer(1, Name, out) ]; net dlnetwork(layers);这里有一个新手容易踩的坑featureInputLayer如果不加 Normalization,none某些版本默认会对输入做归一化处理这会直接改变Fourier特征的频率尺度导致后面怎么调σ都没效果。所以这个参数一定要显式写出来。3.2 采样点生成与Fourier特征映射内部配点在整个时空区域内用均匀随机采样初始条件点固定t0边界点分别固定在x-1和x1。把这些普通数组转成dlarray之后才能参与自动微分。Fourier映射函数的核心就是按坐标拼接后算正弦余弦。%% 2. 采样 % 内部配点x 在 [-1,1]t 在 [0,1] x_f rand(1, Nf) * 2 - 1; t_f rand(1, Nf); % 初始条件点 x_ic rand(1, Nic) * 2 - 1; t_ic zeros(1, Nic); % 边界点 x_bc_l -ones(1, Nbc); t_bc_l rand(1, Nbc); x_bc_r ones(1, Nbc); t_bc_r rand(1, Nbc); % 全部转为 dlarray格式 CB dlXf dlarray(x_f, CB); dlTf dlarray(t_f, CB); dlXic dlarray(x_ic, CB); dlTic dlarray(t_ic, CB); dlXbcL dlarray(x_bc_l,CB); dlTbcL dlarray(t_bc_l,CB); dlXbcR dlarray(x_bc_r,CB); dlTbcR dlarray(t_bc_r,CB); % Fourier 特征映射函数 function z fourierFeature(B, x, t) xt cat(1, x, t); % 2xN z [cos(B * xt); sin(B * xt)]; % (2*m_fourier)xN end这里的技巧是x和t在传给模型之前不拼在一起而是在Fourier映射内拼接这样dlarray的梯度能分别流回 x 和 t。如果提前把坐标合并成一个输入后面单独对t求导会麻烦不少。注意MATLAB脚本中局部函数必须放在脚本文件末尾。fourierFeature 和后面要写的 modelLoss 两个 function 统一放到main脚本的最后也可以分别存成 fourierFeature.m 和 modelLoss.m 两个函数文件再在主脚本中直接调用。分块展示只是为了阅读方便别照抄顺序导致编译报错。3.3 自定义损失函数与训练循环modelLoss函数是整套代码的心脏。先用forward拿到配点处的预测值再连续调用dlgradient求 u_t、u_x、u_xx。接着分别计算PDE残差、初始条件误差、边界误差最后用dlgradient(loss, net.Learnables)得到所有网络参数的梯度。%% 3. 自定义损失函数 function [loss, lossPde, lossIc, lossBc, grads] modelLoss(net, B, ... dlXf, dlTf, dlXic, dlTic, dlXbcL, dlTbcL, dlXbcR, dlTbcR, nu) % --- 内部配点上的 PDE 残差 --- uf forward(net, fourierFeature(B, dlXf, dlTf)); ut dlgradient(uf, dlTf); ux dlgradient(uf, dlXf); uxx dlgradient(ux, dlXf); fpde ut uf .* ux - nu * uxx; lossPde mean(fpde .^ 2, all); % --- 初始条件残差 --- uic forward(net, fourierFeature(B, dlXic, dlTic)); eic uic sin(pi * dlXic); % u(x,0) -sin(pi*x) lossIc mean(eic .^ 2, all); % --- 边界条件残差 --- ul forward(net, fourierFeature(B, dlXbcL, dlTbcL)); ur forward(net, fourierFeature(B, dlXbcR, dlTbcR)); lossBc mean(ul .^ 2, all) mean(ur .^ 2, all); loss lossPde lossIc lossBc; grads dlgradient(loss, net.Learnables); end训练循环用adamupdate维护Adam动量每1000轮打印一次各分量损失方便判断哪一项卡住。这里我把每个loss分量都打印出来不是为了好看而是定位问题的关键如果PDE损失不降通常要考虑σ或配点数如果初始损失不降往往需要加大Nic或提高其权重。%% 4. Adam 训练 trailingAvg []; trailingAvgSq []; figure; for iter 1:nIter [loss, lossPde, lossIc, lossBc, grads] dlfeval(modelLoss, net, B, ... dlXf, dlTf, dlXic, dlTic, dlXbcL, dlTbcL, dlXbcR, dlTbcR, nu); [net, trailingAvg, trailingAvgSq] adamupdate(net, grads, ... trailingAvg, trailingAvgSq, iter, lr); if mod(iter, 1000) 0 fprintf(iter%4d loss%.3e pde%.3e ic%.3e bc%.3e\n, ... iter, extractdata(loss), extractdata(lossPde), ... extractdata(lossIc), extractdata(lossBc)); semilogy(iter, extractdata(loss), b.); hold on; drawnow; end end这里要特别提醒一个小细节在modelLoss内部用forward做前向计算而不是predict。因为在自定义训练循环里predict在某些版本下不会保留完整计算图dlgradient会拿不到梯度。这个坑我周围至少两个人踩过损失函数里报错“Unable to compute gradient”的时候先检查是不是这里写成了predict。3.4 预测与结果可视化训练完网络已经是一个可微函数。在网格上把坐标点铺开批量送入网络预测再把输出reshape成网格形状。用surf画三维曲面可以直观看到初始正弦波如何逐渐演变成激波。%% 5. 结果可视化 [xg, tg] meshgrid(linspace(-1, 1, 200), linspace(0, 1, 100)); xgVec reshape(xg, 1, []); tgVec reshape(tg, 1, []); uPred extractdata(forward(net, fourierFeature(B, ... dlarray(xgVec, CB), dlarray(tgVec, CB)))); uPred reshape(uPred, size(xg)); figure; surf(xg, tg, uPred, EdgeColor, none); xlabel(x); ylabel(t); zlabel(u); title(PINN Fourier Feature 预测的 Burgers 方程解);到这里一套完整可运行的代码就齐了。整体流程是生成坐标 → Fourier特征映射 → 网络前向 → 自动微分算残差 → 三项损失求和 → 反向传播更新参数。把这个框架理解透之后再换成其他方程也就是改残差表达式的事。4. 调参经验与常见问题排查代码能跑通只是第一步实际训练中一定会碰到各种状况。下面把我在调这个算例时踩过和见过的坑集中整理一下。4.1 训练出问题的排查顺序按照“先看损失分量再调对应参数”的思路我建议用下面这个速查表快速定位。现象优先排查调整手段总loss不降或下降极慢学习率是否合适、B矩阵σ是否过小尝试lr1e-2或1e-4加大σ训练中出现NaN学习率过大、σ过大、梯度爆炸降低lr、降低σ、减少m_fourier激波区域过于平滑Fourier特征频率不足增大σ、增大m_fourier初始条件附近误差大初始点太少或权重不足增加Nic、在损失中提高初始项权重边界处翘起边界点太少增加Nbc并在边界附近加密集点训练后期loss反复波动固定采样点过拟合每500轮重新采样一批配点排查时一定把三个损失分量分别打出来看不要只看总loss。很多情况下PDE残差已经降得很低但初始条件loss还挂着说明网络学会了满足方程却没记住起点这种情况光调内部配点是没用的。4.2 激波附近振荡的具体调法Burgers方程这个算例痛点永远集中在激波处。如果你发现激波位置预测出现了马蹄形振荡或者干脆被抹平了多数是频率带宽不够。第一步把σ从5提高到10第二步把m_fourier从100提高到150两步还不够就检查一下训练有没有收敛到更小损失或者把隐藏层宽度从50扩到80。还有一种非常有效但容易被忽略的手段对激波可能出现的区域做局域加密采样。Burgers方程的激波在时间和空间上都有迹可循t越靠近1陡峭区域越集中在x0附近。在 x∈[-0.3,0.3]、t∈[0.6,1] 这个区域多撒一些内部配点网络会更容易把能量花在刀刃上。不过要记得这算一种“利用了先验知识”的做法如果做研究对比要在论文里如实说明。4.3 MATLAB版本、中文注释与调试细节这个项目依赖Deep Learning Toolbox中的dlnetwork、dlgradient、dlfeval和adamupdate推荐R2021a及以上版本太老的版本连featureInputLayer都没有。如果没有featureInputLayer也可以自己把Fourier特征展开后直接作为网络输入绕开这一层但代码会稍微多几行。关于中文注释乱码这个问题在R2023a/R2023b上特别常见多半是脚本文件编码不是UTF-8。解决办法有两个要么在编辑器里把文件另存为UTF-8要么注释直接写成英文。我自己的习惯是代码注释用英文核心说明写在博客正文里一劳永逸。调试时还有一个很实用的技巧把配点数量先降到2000迭代降到1000跑通全流程确认没有报错再恢复到完整配置。不要一上来就全量训练不然报错时根本分不清是代码问题还是参数问题。每次改动只动一个变量比如只调σ或只调学习率然后把三个损失分量记下来对比这样调参才不是玄学。最后说点我自己的体会。这套代码我前前后后跑过很多次最大的感触是Fourier Feature对PINN的提升不是锦上添花而是决定性的。第一次不加这层时激波附近始终糊成一片加上之后整个解的结构都清晰了这种对比在训练图上体现得非常直接。目前很多研究也开始把B矩阵从固定改成可学习或者更进一步引入贝叶斯PINN做不确定性估计这些都是在固定Fourier Feature基础上的自然延伸。建议你拿代码先跑一遍原配置再动手改σ、改网络宽度把每个损失分量拆开观察收获会比只看结论大得多。
