企业官网建设流程全解析

做数值计算的人，多半会绕不开一维扩散方程。虽然这是最简单的抛物型偏微分方程，但它的求解思路直接决定你后面处理复杂 CFD、传热、地下水溶质运移问题时会不会被数值稳定性折磨。我最早认真做这个题目，是在 MATLAB 里把 FE（前向欧拉）、BE（后向欧拉）和 CN（Crank-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)，其中 (i=0,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{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2} ]

这一格式的截断误差是 (O(\Delta x^2))，也就是说空间离散本身是二阶精度。接下来怎么推进时间，就分出了 FE、BE、CN 三种路线。注意，这里用的是 (n) 时刻的值，还是在 (n+1) 时刻的值，还是两个时刻的平均值，这就是三种格式最本质的区别。

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^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\Delta x^2} ]

整理后得到显式迭代式：

[ u_i^{n+1} = u_i^n + r\bigl(u_{i+1}^n - 2u_i^n + u_{i-1}^n\bigr) ]

所谓“显式”，是因为新的时间层 (n+1) 可以直接用旧层 (n) 的值逐点算出来，不需要求解方程组。写代码最省事，但天下没有免费的午餐。稳定性分析会告诉你，FE 的稳定条件是 (r \le 0.5)。一旦违反，解的振幅会指数放大，数值解直接“爆炸”。这就是为什么 FE 适合用来演示数值稳定性，但在实际工程里用得反而不多。

2.2 BE：后向欧拉格式的隐式求解

后向欧拉（Backward Euler，简称 BE）把空间导数放在 (n+1) 时刻求值：

[ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} ]

整理后，未知量 (u_i^{n+1}) 出现多个，但都在同一层，所以必须联立求解。写成矩阵形式：

[ (I - rA),u^{n+1} = u^n ]

其中 (A) 是空间二阶差分的三对角矩阵。这里的 (I) 是单位矩阵，(u^n) 表示内部节点组成的列向量。BE 的优点是无条件稳定，即使 (r) 再大数值解也不会发散。代价是每个时间步都要解一次线性方程组，而且时间方向只有一阶精度。

第一次写隐式格式的人，最容易被“同一个时间层有多个未知量”搞晕。其实理解起来很简单：你把未知量放在等式左边，已知量放在等式右边，剩下的事情就是解一个线性系统。在 MATLAB 里，这一步通常用\运算完成，不需要自己写高斯消元。对稀疏的三对角系统，这个求解过程非常快。

2.3 CN：Crank-Nicolson 的半隐式平均

Crank-Nicolson（简称 CN）的思路很直接：既然 BE 只用 (n+1) 时刻空间差分，FE 只用 (n) 时刻空间差分，那取两者的平均是不是更合理？于是有：

[ \frac{u_i^{n+1} - u_i^n}{\Delta t} = \frac{\alpha}{2}\left[ \frac{u_{i+1}^{n} - 2u_i^{n} + u_{i-1}^{n}}{\Delta x^2} + \frac{u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}}{\Delta x^2} \right] ]

整理成矩阵形式后变成：

[ \left(I - \frac{r}{2}A\right)u^{n+1} = \left(I + \frac{r}{2}A\right)u^n ]

CN 的好处是时间方向达到二阶精度，且无条件稳定。数值实践里，它通常是最平衡的选择：精度高一点，稳定性也不用提心吊胆。代价是同样需要解线性方程，而且矩阵构造比 BE 稍微复杂一点点。

很多人问：CN 是不是就是把 BE 和 FE 做算术平均？从形式上看确实有这种直觉，但从截断误差的角度看，CN 是把时间中点 (t_{n+1/2}) 处的空间导数做中心近似，所以时间方向直接变成二阶精度。这个“半显半隐”的构造，让它同时拿到了稳定性和精度的好处。

2.4 三者在精度和稳定性上的直接对比

用一张表总结会比较直白：

这里面有个很容易踩的误区：无条件稳定不等于结果一定准确。BE 和 CN 虽然不“爆炸”，但如果 (\Delta t) 取得过大，时间离散误差会严重失真，甚至出现非物理振荡。稳定性只是说误差不会被放大到无穷，不保证误差小。这就像开车不超速不代表一定安全，路况不好照样可能翻车。

3. Matlab 实现全过程

3.1 参数设置与初始条件

实现的第一步，是把模型参数、网格参数和初始条件准备好。以区间 ([0,1])、扩散系数 (\alpha=0.02)、总时间 (T=0.5) 为例，典型的 MATLAB 设置如下：

L = 1; % 空间长度 T = 0.5; % 总时间 N = 100; % 空间网格数 M = 500; % 时间步数 alpha = 0.02; % 扩散系数 dx = L / N; x = linspace(0, L, N+1)'; % 所有节点坐标 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)，而不是 (N+1)。因为两端边界值已知，不需要进未知向量。如果把边界节点也算进去，矩阵会变成不可直接求解的形式，或者需要额外处理边界行。

接着基于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: 右侧矩阵

命名很关键：对 FE，A_FE是乘在旧时间层上的；对 BE 和 CN，左侧矩阵是用来解线性方程组的系数矩阵，右侧矩阵是乘在旧时间层上的。建议命名时把left、right写清楚，否则隔两天再看代码，你自己都会分不清哪个是哪个。

3.3 时间推进循环与结果存储

时间循环本身并不复杂，核心是每个时间步按格式更新内部节点向量u_inner：

method = 'CN'; % 可改为 'FE' 或 'BE' u_inner = IC(2:N); % 只取内部节点初始值 U_store = zeros(N+1, M+1); 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, n+1) = u; % 内部节点存起来 end

注意method要在循环之前定义成字符串变量。FE 不需要解线性方程组，是一个稀疏矩阵乘法；BE 和 CN 则用\运算求解稀疏线性系统。对几千个节点的规模，\的耗时完全在可接受范围内。

边界节点之所以不用赋值，是因为U_store(1,:)和U_store(N+1,:)一直保持 0。如果边界值非零，只需要在循环后对首尾两行赋值即可。另一种更稳妥的做法是，在每一步循环后都显式更新首尾节点：

U_store(1, n+1) = 0; % 左边界 u(0,t) U_store(N+1, n+1) = 0; % 右边界 u(L,t)

这样即使你已经初始化了全部节点，也不会因为某次操作意外覆盖边界产生错误。

3.4 可视化与误差分析

数值解算好后，最好直接画三维曲面图，可以非常直观地看到扩散过程中曲线逐渐被拉平的过程：

t = linspace(0, T, M+1); 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 标准算例设计

为了让三种方法在同一起跑线比较，我建议固定空间网格，只改变时间步数，观察时间方向的收敛行为。以 (N=100)、(\alpha=0.02)、(T=0.5) 为例，取 (M=200,400,800,1600)，对应的 (r) 是：

[ r = \frac{0.02 \times 0.5/M}{(1/100)^2} = \frac{100}{M} ]

所以 (M=200) 时 (r=0.5)，(M=400) 时 (r=0.25)，依次减半。这样做的好处是三个方法都在稳定范围内，可以公平比较时间步长对误差的影响。

这里要注意：如果直接取 (M=50)，(r=2)，FE 会当场爆炸，得到 (10^{30}) 量级的数值，根本无法放进同一张表。做对比实验时，要么先把 FE 的稳定性限制考虑进去，要么明确告诉大家“本组实验只用于展示显式格式的失稳”。

4.2 收敛阶验证

从数值结果看，FE 和 BE 的全局时间误差随 (\Delta t) 减小大体是线性下降，也就是斜率接近 1；CN 的误差下降明显更快，斜率接近 2。如果你把误差取对数后做最小二乘拟合，得到的直线斜率就是数值实验中的收敛阶。

下面是一组示意数据，展示了三种方法在 (M=200,400,800,1600) 下的最大误差变化趋势：

这个表不是精确数据，只是用来示意趋势。关键是你会看到 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，其他参数不变，此时 (r=0.5)，FE 勉强稳定。再把 (M) 改成 80，(r=1.25)，重新运行 FE，你会看到数值解在几步之内就出现高频振荡，再过几步直接变成 (10^{30}) 量级的巨大数值。这就是显式条件稳定的现场教学。

比较有意思的是，同样的参数下，BE 和 CN 还能正常算下去。你虽然能算，但要注意解可能被“磨平”了：BE 对高频分量衰减特别厉害，当时间步较大时，解会变得过于平滑，峰值提前消失。CN 相对好很多，但也会在初始阶段出现轻微振荡，尤其是初始条件含有尖锐变化时。

如果想让爆炸过程看得更清楚，可以在循环里实时记录每一步的峰值：

peak(n) = max(abs(u));

当峰值开始指数增长时，说明格式已经不稳定。这是排查代码问题的有力工具。

5. 代码调试和常见问题

5.1 矩阵尺寸对不上怎么办

这是新手最容易遇到的问题。空间节点数是 (N+1)，内部节点数是 (N-1)。矩阵 (A) 的维度必须是 (N-1)，所以如果代码里写成A = spdiags(... , N+1, N+1)，再拿它去乘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 个内部节点的方程是：

[ (1+2r)u_1^{n+1} - r u_2^{n+1} = 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 往往性价比最高，因为它同时改善了时间精度和无条件稳定性。

一个经验公式是：对于瞬态传热模拟，先取 (r=0.5) 算一次，然后把 (r) 减半再算一次，两次结果差异小于 1% 就认为网格和时间步长足够。如果差异很大，继续减小 (r) 直到收敛。

5.4 为什么数值解在初期有振荡

即使 CN 是二阶精度，当初始条件不光滑或者 (r) 较大时，解的傅里叶分量中高频成分在一开始会被放大或衰减不充分，形成小幅振荡。这个现象常见于方波初始条件或者分段常数初值。处理办法通常是把 (\Delta t) 调小，或者用带限制器的方法。对简单扩散问题，我一般优先把 (r) 控制在 1 以内，振荡基本就不明显了。

还有一种情况是边界突变导致的振荡。比如初始条件在边界处不为零，但 Dirichlet 边界强制它为零，这种不连续会造成初期剧烈变化。解决办法是让初始条件在边界附近光滑过渡，或者用更小的初始时间步。

5.5 如何确认代码没有隐藏bug

我曾经被一个“看似正确、结果却完全不对”的代码折磨过很久。排查技巧主要有三个：第一，用解析解做对比，这是最可靠的；第二，检查能量守恒或无界增长趋势，扩散方程的解不应随时间增大；第三，做网格收敛性测试，如果加密网格后结果明显变化，说明当前网格太粗。

另外，代码写完后可以先用一个极小的例子手算一两步。比如 (N=4)、(M=2)，用笔算验证第一个时间步的结果，再和 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，尺寸为 (N+1)*(M+1) % 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) 值、误差、和版本号记在一起，方便回头复现。别小看这个习惯，等过两个月再翻代码时，你会感谢当时的自己。