做数值计算的人多半会绕不开一维扩散方程。虽然这是最简单的抛物型偏微分方程但它的求解思路直接决定你后面处理复杂 CFD、传热、地下水溶质运移问题时会不会被数值稳定性折磨。我最早认真做这个题目是在 MATLAB 里把 FE前向欧拉、BE后向欧拉和 CNCrank-Nicolson三种有限差分格式一一写出来然后拿同一组初始条件去对比精度和稳定性。说实话公式推导一个下午就懂了但真正在 MATLAB 里把三对角矩阵、时间循环和边界条件组装起来反而花了我好几天。这篇文章就把这套 1D 扩散方程的 Matlab 实现完整拆开讲清楚从方程离散、格式推导到代码落地、结果对比和调试经验一次性给你说明白。适合正在学数值方法的本科生、需要快速搭扩散模型的研究生以及所有想用有限差分求解抛物型方程的初学者。1. 问题设定与有限差分基础1.1 一维扩散方程到底在描述什么一维扩散方程的标准形式是[ \frac{\partial u}{\partial t} \alpha \frac{\partial^2 u}{\partial x^2}, \quad 0 x L, ; t 0 ]其中 (u(x,t)) 是待求的物理量(\alpha) 是扩散系数可以理解为热传导系数或分子扩散系数。这个方程描述的是物理量在浓度梯度、温度梯度等作用下从一个高值区域慢慢向低值区域传递最终被“拉平”的过程。初始条件 (u(x,0)f(x)) 给出初始分布边界条件则说明两端是固定值、固定通量还是辐射边界。为什么一维方程值得单独拿出来说因为它是理解抛物型方程数值方法的“最小模型”。把空间离散和时间推进的逻辑在一维搞清楚二维、三维、变系数、非线性项大部分都是在这个框架上做加法。如果一维扩散方程的有限差分代码你调不明白后面学 ADI、对流扩散、反应扩散之类的东西会非常吃力。反过来只要这个基础打牢很多复杂问题都能举一反三。1.2 空间网格和时间步进的基本做法有限差分法的核心是把连续偏导数替换成离散差商。先把空间区间 ([0,L]) 剖成 (N) 段每段长度 (\Delta x L/N)节点坐标 (x_i i\Delta x)其中 (i0,1,\dots,N)。时间区间 ([0,T]) 剖成 (M) 步时间步长 (\Delta t T/M)记 (t_n n\Delta t)。数值解记为 (u_i^n \approx u(x_i,t_n))。空间二阶导数用中心差分近似[ \frac{\partial^2 u}{\partial x^2}\bigg|{x_i} \approx \frac{u{i1}^n - 2u_i^n u_{i-1}^n}{\Delta x^2} ]这一格式的截断误差是 (O(\Delta x^2))也就是说空间离散本身是二阶精度。接下来怎么推进时间就分出了 FE、BE、CN 三种路线。注意这里用的是 (n) 时刻的值还是在 (n1) 时刻的值还是两个时刻的平均值这就是三种格式最本质的区别。1.3 先记住一个关键参数 r在写代码之前一定要先把公式整理成无量纲时间步参数。记[ r \frac{\alpha \Delta t}{\Delta x^2} ]这个参数在很多教材里叫网格傅里叶数也有人直接叫扩散数。它的物理含义是在一个时间步里扩散信息能跨越多少个空间网格。如果 (r) 太大意味着每个时间步信息跨过的网格太多显式格式很快就会出问题。后面所有稳定性和精度讨论最终都会落到 (r) 怎么选上。这个参数的重要性可以类比成“信号传播速度”数值方法总是希望在一个时间步内信息不要移动太多网格否则就容易失真。对 FE 来说(r) 超过 0.5 就直接完蛋对 BE 和 CN 虽然没有硬性上限但 (r) 太大会让时间离散误差变大解失真甚至振荡。所以不管用什么格式第一步先算 (r)养成这个习惯能少踩很多坑。2. 三种时间格式的原理与区别2.1 FE前向欧拉格式的显式推进前向欧拉Forward Euler简称 FE是思路最简单的一种时间导数用前向差分空间导数用 (n) 时刻的值代入。离散方程写成[ \frac{u_i^{n1} - u_i^n}{\Delta t} \alpha \frac{u_{i1}^n - 2u_i^n u_{i-1}^n}{\Delta x^2} ]整理后得到显式迭代式[ u_i^{n1} u_i^n r\bigl(u_{i1}^n - 2u_i^n u_{i-1}^n\bigr) ]所谓“显式”是因为新的时间层 (n1) 可以直接用旧层 (n) 的值逐点算出来不需要求解方程组。写代码最省事但天下没有免费的午餐。稳定性分析会告诉你FE 的稳定条件是 (r \le 0.5)。一旦违反解的振幅会指数放大数值解直接“爆炸”。这就是为什么 FE 适合用来演示数值稳定性但在实际工程里用得反而不多。2.2 BE后向欧拉格式的隐式求解后向欧拉Backward Euler简称 BE把空间导数放在 (n1) 时刻求值[ \frac{u_i^{n1} - u_i^n}{\Delta t} \alpha \frac{u_{i1}^{n1} - 2u_i^{n1} u_{i-1}^{n1}}{\Delta x^2} ]整理后未知量 (u_i^{n1}) 出现多个但都在同一层所以必须联立求解。写成矩阵形式[ (I - rA),u^{n1} u^n ]其中 (A) 是空间二阶差分的三对角矩阵。这里的 (I) 是单位矩阵(u^n) 表示内部节点组成的列向量。BE 的优点是无条件稳定即使 (r) 再大数值解也不会发散。代价是每个时间步都要解一次线性方程组而且时间方向只有一阶精度。第一次写隐式格式的人最容易被“同一个时间层有多个未知量”搞晕。其实理解起来很简单你把未知量放在等式左边已知量放在等式右边剩下的事情就是解一个线性系统。在 MATLAB 里这一步通常用\运算完成不需要自己写高斯消元。对稀疏的三对角系统这个求解过程非常快。2.3 CNCrank-Nicolson 的半隐式平均Crank-Nicolson简称 CN的思路很直接既然 BE 只用 (n1) 时刻空间差分FE 只用 (n) 时刻空间差分那取两者的平均是不是更合理于是有[ \frac{u_i^{n1} - u_i^n}{\Delta t} \frac{\alpha}{2}\left[ \frac{u_{i1}^{n} - 2u_i^{n} u_{i-1}^{n}}{\Delta x^2} \frac{u_{i1}^{n1} - 2u_i^{n1} u_{i-1}^{n1}}{\Delta x^2} \right] ]整理成矩阵形式后变成[ \left(I - \frac{r}{2}A\right)u^{n1} \left(I \frac{r}{2}A\right)u^n ]CN 的好处是时间方向达到二阶精度且无条件稳定。数值实践里它通常是最平衡的选择精度高一点稳定性也不用提心吊胆。代价是同样需要解线性方程而且矩阵构造比 BE 稍微复杂一点点。很多人问CN 是不是就是把 BE 和 FE 做算术平均从形式上看确实有这种直觉但从截断误差的角度看CN 是把时间中点 (t_{n1/2}) 处的空间导数做中心近似所以时间方向直接变成二阶精度。这个“半显半隐”的构造让它同时拿到了稳定性和精度的好处。2.4 三者在精度和稳定性上的直接对比用一张表总结会比较直白格式时间精度空间精度稳定性每个时间步的计算成本FE一阶二阶需要 (r \le 0.5)一次矩阵乘法BE一阶二阶无条件稳定一次线性方程组求解CN二阶二阶无条件稳定一次线性方程组求解这里面有个很容易踩的误区无条件稳定不等于结果一定准确。BE 和 CN 虽然不“爆炸”但如果 (\Delta t) 取得过大时间离散误差会严重失真甚至出现非物理振荡。稳定性只是说误差不会被放大到无穷不保证误差小。这就像开车不超速不代表一定安全路况不好照样可能翻车。3. Matlab 实现全过程3.1 参数设置与初始条件实现的第一步是把模型参数、网格参数和初始条件准备好。以区间 ([0,1])、扩散系数 (\alpha0.02)、总时间 (T0.5) 为例典型的 MATLAB 设置如下L 1; % 空间长度 T 0.5; % 总时间 N 100; % 空间网格数 M 500; % 时间步数 alpha 0.02; % 扩散系数 dx L / N; x linspace(0, L, N1); % 所有节点坐标 dt T / M; r alpha * dt / dx^2; % 网格傅里叶数 fprintf(r %.3f\n, r);初始条件建议选择有解析解的形式方便后续验证误差。最简单的就是正弦函数IC sin(pi * x); % u(x,0) sin(pi*x)如果边界是 Dirichlet 零边界那么两端节点的值始终为 0不需要在求解过程中更新。真正参与计算的是内部节点也就是下标从 2 到 N 的位置。3.2 三对角矩阵组装内部节点数量为 (N-1)。空间二阶差分矩阵 (A) 是一个三对角矩阵对角元素为 (-2)两条次对角元素为 1。MATLAB 里用spdiags构造最合适Nin N - 1; e ones(Nin, 1); A spdiags([e -2*e e], -1:1, Nin, Nin);spdiags的第一个输入是列向量组成的矩阵第二个输入是对应对角线位置。-1:1表示主对角线以及上下两条次对角线。这样得到的A是稀疏矩阵当网格数达到几千、几万时内存占用比diag构造的稠密矩阵小得多线性求解速度也快很多。这里特别提醒A的维度必须是 (N-1)而不是 (N1)。因为两端边界值已知不需要进未知向量。如果把边界节点也算进去矩阵会变成不可直接求解的形式或者需要额外处理边界行。接着基于A构造三种格式的推进矩阵I speye(Nin); A_FE I r * A; % FE: 显式推进矩阵 M_BE I - r * A; % BE: 隐式左侧矩阵 M_CN_left I - 0.5 * r * A; % CN: 左侧矩阵 M_CN_right I 0.5 * r * A; % CN: 右侧矩阵命名很关键对 FEA_FE是乘在旧时间层上的对 BE 和 CN左侧矩阵是用来解线性方程组的系数矩阵右侧矩阵是乘在旧时间层上的。建议命名时把left、right写清楚否则隔两天再看代码你自己都会分不清哪个是哪个。3.3 时间推进循环与结果存储时间循环本身并不复杂核心是每个时间步按格式更新内部节点向量u_innermethod CN; % 可改为 FE 或 BE u_inner IC(2:N); % 只取内部节点初始值 U_store zeros(N1, M1); U_store(:,1) IC; u u_inner; for n 1:M switch method case FE u A_FE * u; case BE u M_BE \ u; case CN u M_CN_left \ (M_CN_right * u); end U_store(2:N, n1) u; % 内部节点存起来 end注意method要在循环之前定义成字符串变量。FE 不需要解线性方程组是一个稀疏矩阵乘法BE 和 CN 则用\运算求解稀疏线性系统。对几千个节点的规模\的耗时完全在可接受范围内。边界节点之所以不用赋值是因为U_store(1,:)和U_store(N1,:)一直保持 0。如果边界值非零只需要在循环后对首尾两行赋值即可。另一种更稳妥的做法是在每一步循环后都显式更新首尾节点U_store(1, n1) 0; % 左边界 u(0,t) U_store(N1, n1) 0; % 右边界 u(L,t)这样即使你已经初始化了全部节点也不会因为某次操作意外覆盖边界产生错误。3.4 可视化与误差分析数值解算好后最好直接画三维曲面图可以非常直观地看到扩散过程中曲线逐渐被拉平的过程t linspace(0, T, M1); mesh(x, t, U_store); xlabel(x); ylabel(t); zlabel(u); colorbar;如果要定量对比三种格式的精度可以用解析解。对 (u(x,0)\sin(\pi x))、零边界条件解析解是[ u_{\text{exact}}(x,t) e^{-\alpha \pi^2 t} \sin(\pi x) ]最后时刻的最大误差或二范数误差可以这样算u_exact exp(-alpha * pi^2 * T) * sin(pi * x); err max(abs(U_store(:,end) - u_exact)); fprintf(最大误差 %.6e\n, err);有了误差值就可以做收敛阶实验了。这里额外提一句画图时如果发现曲面边缘出现锯齿状波纹多半是时间步长过大先不要急着改算法把 (r) 压小再看。4. 数值实验结果对比4.1 标准算例设计为了让三种方法在同一起跑线比较我建议固定空间网格只改变时间步数观察时间方向的收敛行为。以 (N100)、(\alpha0.02)、(T0.5) 为例取 (M200,400,800,1600)对应的 (r) 是[ r \frac{0.02 \times 0.5/M}{(1/100)^2} \frac{100}{M} ]所以 (M200) 时 (r0.5)(M400) 时 (r0.25)依次减半。这样做的好处是三个方法都在稳定范围内可以公平比较时间步长对误差的影响。这里要注意如果直接取 (M50)(r2)FE 会当场爆炸得到 (10^{30}) 量级的数值根本无法放进同一张表。做对比实验时要么先把 FE 的稳定性限制考虑进去要么明确告诉大家“本组实验只用于展示显式格式的失稳”。4.2 收敛阶验证从数值结果看FE 和 BE 的全局时间误差随 (\Delta t) 减小大体是线性下降也就是斜率接近 1CN 的误差下降明显更快斜率接近 2。如果你把误差取对数后做最小二乘拟合得到的直线斜率就是数值实验中的收敛阶。下面是一组示意数据展示了三种方法在 (M200,400,800,1600) 下的最大误差变化趋势M误差FE误差BE误差CN2002.31e-32.35e-36.80e-54001.16e-31.18e-31.70e-58005.80e-45.95e-44.20e-616002.90e-42.98e-41.05e-6这个表不是精确数据只是用来示意趋势。关键是你会看到 CN 的误差比 FE/BE 小一到两个数量级且减半时间步长时CN 误差大约除以 4而 FE/BE 大约除以 2。这就把理论上的收敛阶直观验证了。如果你想让收敛阶计算更严谨可以在不同 (M) 下运行代码并把结果保存到数组里M_list [200 400 800 1600]; err_list zeros(size(M_list)); for k 1:numel(M_list) % 重新运行求解器得到 U_store err_list(k) max(abs(U_store(:,end) - u_exact)); end p polyfit(log(1./M_list), log(err_list), 1); fprintf(收敛阶 ≈ %.2f\n, p(1));对 CN得到的p(1)应该在 2 附近FE 和 BE 则在 1 附近。这就是用数值实验验证理论误差阶的标准做法。4.3 稳定性实验让 FE 当场爆炸稳定性这块一定要亲手做一次。把 (M) 改成 200其他参数不变此时 (r0.5)FE 勉强稳定。再把 (M) 改成 80(r1.25)重新运行 FE你会看到数值解在几步之内就出现高频振荡再过几步直接变成 (10^{30}) 量级的巨大数值。这就是显式条件稳定的现场教学。比较有意思的是同样的参数下BE 和 CN 还能正常算下去。你虽然能算但要注意解可能被“磨平”了BE 对高频分量衰减特别厉害当时间步较大时解会变得过于平滑峰值提前消失。CN 相对好很多但也会在初始阶段出现轻微振荡尤其是初始条件含有尖锐变化时。如果想让爆炸过程看得更清楚可以在循环里实时记录每一步的峰值peak(n) max(abs(u));当峰值开始指数增长时说明格式已经不稳定。这是排查代码问题的有力工具。5. 代码调试和常见问题5.1 矩阵尺寸对不上怎么办这是新手最容易遇到的问题。空间节点数是 (N1)内部节点数是 (N-1)。矩阵 (A) 的维度必须是 (N-1)所以如果代码里写成A spdiags(... , N1, N1)再拿它去乘IC(2:N)要么维度不匹配要么解出来的结果全是错的。我建议在组装矩阵后主动加一行检查assert(isequal(size(A), [Nin Nin]), 矩阵维度错误请检查N的取值);类似的初始向量u0_inner IC(2:N)的长度应该是 (N-1)和矩阵维度完全一致。多用一个节点都会造成\运算报错。如果你遇到Matrix dimensions must agree先别急着改代码把所有数组的size打印一遍十有八九是内部节点数没对齐。这个习惯能帮你省下大量调试时间。5.2 边界条件怎么夹进去上面的实现是 Dirichlet 零边界边界节点值固定为 0不需要进矩阵。如果你要处理非零 Dirichlet 边界比如 (u(0,t)u_L)(u(L,t)u_R)有几种做法。最简单的做法是方程离散时把边界值放到右端项中。以 BE 为例第 1 个内部节点的方程是[ (12r)u_1^{n1} - r u_2^{n1} u_1^n r,u_L ]最后一个内部节点同理要加上 (r,u_R)。如果边界值随时间变化还要在每一步重新更新右端项。建议把边界条件单独写成函数不要散落在主循环里。需要处理 Neumann 边界时事情会稍微复杂。比如 (u_x(0,t)0) 表示左边界绝热通常会用虚拟节点或单侧差分把边界节点消去。这个扩展超出本文范围但你掌握 Dirichlet 之后再学会顺畅很多。5.3 时间步长怎么选才合理实际工程里很少会严格计算 von Neumann 稳定性更多是先用解析解或粗略估计定一个 (r)然后做网格独立性验证。FE 必须满足 (r \le 0.5)BE 和 CN 没有稳定性硬限制但建议初始尝试时 (r) 不要超过 5 或 10否则解虽然不发散时间误差也可能大到离谱。如果发现自己算出的解有明显“台阶状”振荡通常是 (r) 过大。解决办法有三个方向减小 (\Delta t)、减小 (\Delta x)、换用 CN。其中换 CN 往往性价比最高因为它同时改善了时间精度和无条件稳定性。一个经验公式是对于瞬态传热模拟先取 (r0.5) 算一次然后把 (r) 减半再算一次两次结果差异小于 1% 就认为网格和时间步长足够。如果差异很大继续减小 (r) 直到收敛。5.4 为什么数值解在初期有振荡即使 CN 是二阶精度当初始条件不光滑或者 (r) 较大时解的傅里叶分量中高频成分在一开始会被放大或衰减不充分形成小幅振荡。这个现象常见于方波初始条件或者分段常数初值。处理办法通常是把 (\Delta t) 调小或者用带限制器的方法。对简单扩散问题我一般优先把 (r) 控制在 1 以内振荡基本就不明显了。还有一种情况是边界突变导致的振荡。比如初始条件在边界处不为零但 Dirichlet 边界强制它为零这种不连续会造成初期剧烈变化。解决办法是让初始条件在边界附近光滑过渡或者用更小的初始时间步。5.5 如何确认代码没有隐藏bug我曾经被一个“看似正确、结果却完全不对”的代码折磨过很久。排查技巧主要有三个第一用解析解做对比这是最可靠的第二检查能量守恒或无界增长趋势扩散方程的解不应随时间增大第三做网格收敛性测试如果加密网格后结果明显变化说明当前网格太粗。另外代码写完后可以先用一个极小的例子手算一两步。比如 (N4)、(M2)用笔算验证第一个时间步的结果再和 MATLAB 输出对比。这个方法虽然笨但能迅速定位是公式推导错、矩阵构造错还是边界处理错。6. 写代码时的一些个人习惯6.1 换网格之前先算一遍 r代码敲完第一件事不是急着跑而是手算一遍 (r) 是否落在合理区间。很多“算出来是错的但也不报错”的案例根源就是 (r) 太大或太小。我习惯在fprintf里把 (r)、(\Delta x)、(\Delta t) 三个值在运行时打印出来看一眼心里才有底。fprintf(dx %.4f, dt %.6f, r %.3f\n, dx, dt, r);这一行代码能救命。当别人拿着一张奇怪的图来问我“哪里出问题了”时我第一句话永远是你的 (r) 是多少6.2 把三种格式封装成函数调试阶段可以把三种格式写在一个脚本里方便对比。但一旦要反复换参数、换初始条件、换边界最好还是封装成函数。比如function U solveDiffusion(method, L, T, N, M, alpha, IC) % 返回所有时刻的数值解矩阵 U尺寸为 (N1)*(M1) % method 可选 FE、BE、CN ... end封装之后数据流清晰后续扩展到二维或非线性问题也更容易。你可以在函数内部用switch method分别处理三种格式主脚本只负责传参。再配合一个函数画图整个项目结构就很干净了。6.3 版本兼容性提醒MATLAB 的稀疏矩阵运算接口已经稳定了很多年spdiags、speye、\这些函数在绝大多数版本里都能用。唯一要留意的是不要把method定义成 MATLAB 内置函数名比如不要叫error或format否则会覆盖默认功能。如果所在环境没有 MATLAB用 Octave 跑这段代码也基本兼容只需要把脚本里的中文注释另存为 UTF-8 编码即可。我个人的体会是一维扩散方程的三种格式真正难的不是公式而是把公式翻译成矩阵运算时那些“差一个下标”的问题。只要耐心把边界、内部节点、(N-1) 维度这三件事厘清后面学 ADI 格式、对流扩散方程、非线性反应扩散方程都会顺很多。希望这篇 Matlab 实现笔记能让你少走点弯路。最后再分享一个习惯每次跑完数值实验我会把对应的 (r) 值、误差、和版本号记在一起方便回头复现。别小看这个习惯等过两个月再翻代码时你会感谢当时的自己。