企业官网建设流程全解析

简介：这份PDF面向学习计算流体力学与数值格式的本科生、研究生及Matlab实践者，围绕二维对流方程在迎风C格式、A格式、Lax-Wendroff格式和时间分裂法下的数值求解展开，重点处理∆t/∆x与∆t/∆y取0.5、1、2时的稳定性讨论，并给出u-x-y曲面图与动画制作思路。压缩包仅含1个PDF文件，大小约25.94MB，正文按问题描述、实验原理、差分格式、边界条件、运行结果、结果分析和附录程序组织，便于顺着报告脉络复现计算流程。已有319人学习，适合作为CFD上机实习或课程设计的参考。读者可从中获得四种格式的离散表达式、入流出流边界处理方式、不同Courant数下的结果对比与稳定性分析，以及可直接借鉴的Matlab程序结构和绘图方法，用于核对自身代码与实验现象。

1. 同一个方程、四种写法、三种结局：稳定性不是玄学

二维对流方程 ∂u/∂t + ∂u/∂x + ∂u/∂y = 0，初值是一个中心 8×8 的方形平台，网格 Δx = Δy = 1，计算域 [-16,16]²，推进 90 步。这是计算流体力学里最简单的一个 CFD 算例，Matlab 三十行就能写完。可一旦把 γ = Δt/Δx 分别取 0.5、1、2，四种差分格式给出的结果完全不是一回事：迎风 C 格式在 0.5 稳住、在 1 和 2 处爆掉；A 格式三档全爆；Lax-Wendroff 也全爆，要把 γ 压到 0.3 才看得见一个像样的波形；时间分裂法下的 C 格式却在 0.5 和 1 都收敛，直到 2 才散。

四个格式写的是同一个方程，代码骨架几乎一样，差别只在模板取哪几个邻居、系数乘多少。所以「稳定性」从来不是调参调出来的玄学，而是格式本身写死了的一条硬边界。搞懂这条边界怎么来的，才算真正把 CFD 入门，而不是把有限差分当黑匣子调着玩。

2. 放大因子 G：冯·诺依曼分析如何把稳定性压成一道复平面几何题

2.1 误差和真解走同一个差分方程，这是整套推导的地基

差分格式可以写成 u^{n+1} = L(u^n) 的形式，L 是一个线性算子。设数值解为 U，精确解为 u，误差 ε = U - u。因为方程和格式都是线性的，且两者满足同样的边界条件，把两个式子相减，ε 满足的差分方程和 U 满足的完全一致——只是初值从平台函数换成了初始误差分布。

这一步是整个稳定性分析的基础：误差不受解答本身影响，它会独立地按格式规定的规律演化。既然误差自己会走，那么判断一个格式行不行，就等价于判断「任意初始误差在反复施加 L 之后会不会被放大到无穷」。如果存在某个误差模态每次迭代都被放大一点点，90 步、900 步之后就会盖过真解，我们在图上看到的就是数值爆炸而不是波形传播。

常见的实现思路是只在常系数、周期边界下讨论这件事，因为这时 L 是一个循环矩阵，特征函数是离散傅里叶模态，误差可以干净地分解成一堆平面波的叠加。这也是冯·诺依曼分析的适用前提，绝大多数 CFD 格式的第一轮筛选都靠它。

2.2 模态代入：把逐点递推换成逐波数缩放

取单个误差模态 ε^n_{i,j} = A_n · exp(i(k1·x_i + k2·y_j))，其中 k1、k2 是两个方向的波数，A_n 是第 n 步的复振幅。把它代进误差传播方程，左右两边 exp 因子会整体约掉，剩下的只有振幅之间的比例关系，定义这个比例为放大因子 G = A_{n+1} / A_n。

以迎风 C 格式为例，误差方程是

ε^{n+1}{i,j} = ε^n{i,j} - γ(ε^n_{i,j} - ε^n_{i-1,j}) - γ(ε^n_{i,j} - ε^n_{i,j-1})

代入后，ε^n_{i-1,j} = ε^n_{i,j}·exp(-i·k1·Δx)，ε^n_{i,j-1} = ε^n_{i,j}·exp(-i·k2·Δy)。记 β = k1·Δx、θ = k2·Δy，整理得到

G = 1 - 2γ + γ·exp(-iβ) + γ·exp(-iθ)

β 和 θ 各自扫过 [-π, π]，就覆盖了网格能分辨的全部波数。稳定性要求对每一个 (β, θ)，|G| ≤ 1 都成立；只要存在一组波数让 |G| > 1，对应模态就会被指数放大，格式不稳定。这个结论在常系数周期边界下是充要的，作为工程判据足够用。

用几行 Matlab 可以把某个格式的 |G| 极值直接扫出来，比手推快得多，也方便验证推导有没有算错：

% 扫描全部波数，求放大因子的最大模（迎风 C 格式 / A 格式 / LW 格式通用） gamma = 0.5; % gamma = dt/dx = dt/dy beta = linspace(-pi, pi, 601); theta = linspace(-pi, pi, 601); [B, T] = meshgrid(beta, theta); G_C = 1 - 2*gamma + gamma*exp(-1i*B) + gamma*exp(-1i*T); % 迎风 C G_A = 1 - 1i*gamma*(sin(B) + sin(T)); % A 格式 G_LW = G_A - 2*gamma^2*(sin(B/2).^2 + sin(T/2).^2); % Lax-Wendroff fprintf('C : %.6f\n', max(abs(G_C (:)))); fprintf('A : %.6f\n', max(abs(G_A (:)))); fprintf('LW : %.6f\n', max(abs(G_LW(:))));

gamma 就是网格雷诺数意义上的库朗数，B、T 是两个方向的无量纲波数。把 gamma 从 0.1 扫到 2，打印三条曲线，你会看到 C 格式的 |G|max 在 γ = 0.5 处从 1 拐头向上，A 格式从 γ 稍微大于 0 开始就永远大于 1，LW 格式的拐点在 0.3536 附近。三条曲线，三种命运，后面几章就是把它们各自的原因拆开。

2.3 |G|≤1 的代数与几何双重读法

把迎风 C 格式的 G 写成矢量叠加的形式最直观：

G = (1 - 2γ, 0) + (γ·cosβ, -γ·sinβ) + (γ·cosθ, -γ·sinθ)

第一个矢量是实轴上长度 1 - 2γ 的线段，后两个是模长都为 γ、角度分别由 β 和 θ 决定的旋转矢量。当 β = θ 时，两个旋转矢量同向，它们的水平分量叠加、垂直分量也叠加，这是最容易冲出去的方向：两个 γ 长的箭头并列甩向左边，正好抵消掉实轴段的收缩，也最容易越过单位圆。

也可以老老实实算一遍。令 β = θ、c = cosβ、s = sinβ：

|G|² = (1 - 2γ + 2γc)² + (2γs)² = 1 - 4γ + 8γ² + 4γ(1 - 2γ)·c

γ ≤ 0.5 时系数 1 - 2γ ≥ 0，最大值在 c = 1 处取到，代入得 |G|²max = 1，恰好压在单位圆上；γ > 0.5 时系数变负，最大值挪到 c = -1，|G|²max = (1 - 4γ)²，模长 |G|max = |1 - 4γ|。γ = 1 时是 3，γ = 2 时是 7。这意味着每推进一个时间步，最快那个误差模态就放大 7 倍，90 步之后是 7⁹⁰ 量级，surf 图上只剩几根戳破坐标轴的尖刺，其余部分被压成一张平板。

注意：发散在图上的表现形式常常不是「数值变成 NaN」，而是「极少数点数值巨大、其余点看起来还是零」。看到这种形态先别怀疑初值，优先怀疑放大因子越界。

3. 迎风 C 格式：γ≤0.5 这条线究竟卡在哪个方向

3.1 依赖域：迎风模板把信息流向写进了网格

迎风 C 格式的完整形式是

u^{n+1}{i,j} = u^n{i,j} - γ(u^n_{i,j} - u^n_{i-1,j}) - γ(u^n_{i,j} - u^n_{i,j-1})

它只取 (i,j)、(i-1,j)、(i,j-1) 三个点。这三个点正好位于 (i,j) 点的上游：特征线从西北方向来，信息沿 x、y 正方向传播，模板就往负方向偏一格。这种「偏」不是随意选择，而是把 PDE 的特征方向编码进离散模板，让数值依赖域尽量贴近真实依赖域。

一维情况下，CFL 条件要求特征线在 Δt 内跨过的距离不超过一个网格，即 γ ≤ 1。二维时两个方向同时推进，Δt 内波在 x 方向走 γ 格、y 方向走 γ 格，加起来不能超过一个网格的对角尺度，于是条件变成 γ_x + γ_y ≤ 1。取 γ_x = γ_y = γ，就得到 γ ≤ 0.5。这跟上一章代数推导的结果完全吻合——两条路，同一个答案。

迎风 C 格式是一阶精度的，代价是数值耗散明显：课件里那个方形平台在 γ = 0.5 下推进 90 步后，边缘被抹圆、峰值被削低，看起来像被风吹散的一团雾。但耗散的另一面是稳：它老老实实把高频误差压下去，不会反过来放大。

3.2 入流给值、出流外插：边界处理决定了波形能不能走出去

边界条件是这类算例里最容易翻车的地方，比格式本身还容易。波向 x、y 正方向传播，那么左边界（i = 1）和下边界（j = 1）是入流边界，物理上信息从外面流进来，直接给定值 u = 0；右边界（i = end）和上边界（j = end）是出流边界，信息从域内流出去，不能强行给定值，否则波撞上墙会反弹回来，形成一圈假的回波。

常见的处理是线性外插，把内部最后两点的斜率延长一格：

u_{i+1,j} = 2u_{i,j} - u_{i-1,j} u_{i,j+1} = 2u_{i,j} - u_{i,j-1}

这个式子的含义是「边界外那一格的值，按线性趋势推」——等价于二阶的单侧差分，允许波平滑地离开计算域而不产生明显反射。它的误差量级比格式本身低一阶，不会污染内部解，是出流边界最省事也最常用的做法。

实现顺序要特别注意：内点循环跑完之后再统一补边界，而且近边界那一圈内点的循环范围要相应缩小。如果一边循环一边改边界，后面的点会用到已经更新过的邻居，等价于在格式里混进了半隐式成分，稳定性分析就不再适用。

3.3 Matlab 实现与三档 γ 的对照

下面是迎风 C 格式的主体循环。网格用 meshgrid(x, y) 生成，注意此时 X(i,j) = x(j)、Y(i,j) = y(i)，第一维对应 y、第二维对应 x，方台的索引区间要和这个约定对齐：

dx = 1; dy = 1; gamma = 0.5; % gamma = dt/dx = dt/dy x = -16:dx:16; y = -16:dy:16; [X, Y] = meshgrid(x, y); % X(i,j)=x(j), Y(i,j)=y(i) U1 = zeros(size(X)); U1(13:21, 13:21) = 1; % 初值：x,y 均落在 [-4,4] 的方形平台 U2 = zeros(size(U1)); for n = 1:90 for i = 2:size(X,1) for j = 2:size(X,2) % 迎风 C 格式：只用上游的 (i-1,j) 与 (i,j-1) U2(i,j) = U1(i,j) - gamma*(U1(i,j) - U1(i-1,j)) ... - gamma*(U1(i,j) - U1(i,j-1)); end end U2(1,:) = 0; U2(:,1) = 0; % 入流边界 U2(end,:) = 2*U2(end-1,:) - U2(end-2,:); % 出流：x 向外插 U2(:,end) = 2*U2(:,end-1) - U2(:,end-2); % 出流：y 向外插 U1 = U2; end

gamma 是库朗数，取 0.5 时正好踩在稳定边界上；dx、dy 固定为 1 时它就等于时间步长。索引 13:21 对应 x 从 -16 起算第 13 到第 21 个点，即 [-4, 4] 共 9 个网格点。入流边界的赋值放在内点循环之后，是为了保证边界值覆盖掉循环里可能算出的无效值；出流外插也在同一位置，顺序是先 x 后 y，两者互不依赖。

把 gamma 改成 1 或 2，同一个脚本跑出来的结果就是上一章算过的 3ⁿ 和 7ⁿ。这不是步长没调好，是格式在 γ > 0.5 时对某些波数的 |G| 已经大于 1，任何初值、任何边界处理都救不回来。

4. A 格式、Lax-Wendroff 与时间分裂：γ 从 0 到 1 的稳定性谱系

4.1 A 格式：纯虚特征值撞上显式欧拉稳定域的边界

A 格式把迎风改成中心差分：

u^{n+1}{i,j} = u^n{i,j} - (γ/2)(u^n_{i+1,j} - u^n_{i-1,j}) - (γ/2)(u^n_{i,j+1} - u^n_{i,j-1})

代入模态后，两个 exp 相减变成正弦：

G = 1 - i·γ·(sinβ + sinθ)

实部恒为 1，虚部是一个纯虚数，于是 |G|² = 1 + γ²(sinβ + sinθ)²。除了 sinβ + sinθ = 0 的少数模态，其余模态的模长全部大于 1；γ 越大、偏离越多。γ = 0.5、β = θ = π/2 时，|G|² = 1 + 4×0.25 = 2，单步放大 1.414 倍，90 步后放大约 10¹¹ 量级。

这个结果可以直接在复平面上读懂。把格式写成 G = 1 + z 的形式，z = -iγ(sinβ + sinθ) 落在虚轴上。显式欧拉的绝对稳定域是以 -1 为圆心、半径为 1 的圆盘，它只在原点处与虚轴相切。任何纯虚的 z 都会把 1 + z 推出单位圆，只有一个点除外。换言之，中心差分 + 显式时间推进这个组合，天生就把特征值放在了稳定域之外，跟 γ 取多少无关。A 格式的不稳定是无条件的，不是「参数没调好」。

实践里看到 A 格式在 γ 很小时前几步还算正常，那是初值本身光滑、高频误差还没被激发出来的假象。跑够步数，误差一定会长出来。

4.2 Lax-Wendroff：二阶修正项换来精度，也换来 γ≤0.3536

Lax-Wendroff 的思路是：既然中心差分是一阶时间精度，那就按泰勒展开补一个二阶时间项，用空间二阶导数替换掉它，得到

u^{n+1}{i,j} = u^n{i,j} - (γ/2)(u^n_{i+1,j} - u^n_{i-1,j}) - (γ/2)(u^n_{i,j+1} - u^n_{i,j-1}) + (γ²/2)(u^n_{i+1,j} - 2u^n_{i,j} + u^n_{i-1,j}) + (γ²/2)(u^n_{i,j+1} - 2u^n_{i,j} + u^n_{i,j-1})

前两项和 A 格式一模一样，后两项是新增的阻尼项，形式上就是两个方向的离散拉普拉斯算子。放大因子相应变成

G = 1 - i·γ(sinβ + sinθ) - 2γ²(sin²(β/2) + sin²(θ/2))

虚部由 -iγ(sinβ + sinθ) 贡献，随 γ 线性增长；阻尼项由 -2γ²(...) 贡献，随 γ² 增长。两者在 γ 较小时后半部分压得住前半部分，γ 一大会反过来。平衡点在哪，就决定了稳定域的上限。

二维的情况比一维严格得多。一维 Lax-Wendroff 的稳定条件是 γ ≤ 1，但二维方形网格上，对角波 β = θ 会让两个方向的虚部同相叠加，达到单个方向的两倍，而阻尼项的叠加跟不上这个速度。按已知的判据，二维稳定性条件写成 K1^{2/3} + K2^{2/3} ≤ 1，K1 = K2 = γ 时解出 γ ≤ 2^{-3/2} ≈ 0.3536。这个数比 0.5 还小，比一维的 1 小了近三倍。

数值实验能验证这一点：γ 取 0.5、1、2 三档全炸，把 γ 降到 0.3，波形立刻恢复正常。0.3 < 0.3536，刚好在安全区内。代价也很清楚——同样推进 90 个时间步，γ = 0.3 意味着物理时间只走了 27，而迎风 C 格式在 γ = 0.5 时走了 45。想推进到同样的物理时刻，Lax-Wendroff 要多花将近一倍的时间步，这是二阶精度在二维买的单。

4.3 时间分裂法：两个一维迎风串联，稳定域反而放宽到 γ≤1

时间分裂法换了个思路：不直接离散二维方程，而是把它拆成两个一维方程，先沿 y 推进一次，再用中间结果沿 x 推进一次：

u^{n+1/2}{i,j} = u^n{i,j} - γ(u^n_{i,j} - u^n_{i,j-1}) u^{n+1}{i,j} = u^{n+1/2}{i,j} - γ(u^{n+1/2}{i,j} - u^{n+1/2}{i-1,j})

每一步都是标准的一维迎风格式。因为两步是算子的串联，放大因子就是两个一维放大因子的乘积：

G = [1 - γ(1 - e^{-iθ})] · [1 - γ(1 - e^{-iβ})]

分别算每个因子的模。以 y 方向那个为例，展开得 |·|² = 1 - 2γ + 2γ² + 2γ(1 - γ)cosθ。cosθ = 1 时值恒为 1；cosθ = -1 时值为 (1 - 2γ)²，只要 γ ≤ 1 就不超过 1。中间取值更小。所以单个因子在 0 ≤ γ ≤ 1 内满足 |G| ≤ 1，乘积自然也不超过 1。稳定域就这样从 0.5 放宽到了 1。

原因不复杂：二维迎风 C 格式里存在交叉项，x、y 的推进同时发生，误差在两个方向上的贡献会叠加。分裂法把它们错开了，任何时刻只在一个方向上推进，交叉项根本不存在，只剩两个独立的一维问题。一维迎风的条件是 γ ≤ 1，乘积的条件也是 γ ≤ 1。这就是「拆开反而更稳」的全部秘密。

实现上有一个细节。既然每步只沿一个方向推进，另一步的推进方向在时间上是分开的，那么第一步求出的中间量在第二步里会被完整消费掉，边界处理也要分两次做。代码如下：

r = 1.0; % 分步推进时的系数，取 1 仍稳定 for n = 1:90 % 第一步：沿 y 方向推进（对应 i 方向索引） for i = 2:size(X,1) for j = 1:size(X,2) U3(i,j) = U1(i,j) - r*(U1(i,j) - U1(i-1,j)); end end U3(end,:) = 2*U3(end-1,:) - U3(end-2,:); % 中间步出流外插 % 第二步：沿 x 方向推进（对应 j 方向索引） for i = 1:size(X,1) for j = 2:size(X,2) U2(i,j) = U3(i,j) - r*(U3(i,j) - U3(i,j-1)); end end U2(:,end) = 2*U2(:,end-1) - U2(:,end-2); % 整步出流外插 U2(1,:) = 0; U2(:,1) = 0; % 入流边界 U1 = U2; end

r 是分步推进的系数，理论上每一步用半时间步更严谨，取全步时稳定性判据会相应收紧；算例里直接取 r = 0.5 或 1 都能出结果，这正是分裂法宽容度的体现。U3 是中间层，注意它只在两步之间存活，不要复用成下一轮的初值。两次出流外插必须分开做：第一步的外插作用于 y 方向的边界，第二步的外插作用于 x 方向，混在一起会让两个方向的外插互相污染。

分裂法的代价是时间精度。两个一维算子串联，交换性误差是一个 O(Δt) 量级的分裂误差，因此整体时间精度只有一阶，达不到 Lax-Wendroff 的二阶。稳、准、快三者在这里第一次出现明显的取舍。

4.4 奇偶步交换：为什么要打乱固定的 x-y 顺序

固定「先 y 后 x」的分裂顺序会带来一个隐蔽的问题：边界外插、入流赋值都是按方向做的，两个方向的处理并不对称。长期跑下来，误差会沿对角线方向单向积累，波形出现定向的偏斜——看起来像是格式不稳，其实是分裂顺序带来的各向异性。

解决办法是奇偶步交替：奇数步用 y → x，偶数步用 x → y。这样两个方向的处理在时间上被平均掉，一阶的分裂误差不会持续往同一个方向堆。算例里明确写了这个安排，代码结构就是把两个分支并排放在 if mod(n,2) 的两侧，各自带一套完整的边界处理。

if mod(n,2) == 1 % 奇数步：先 y 后 x % ... y 方向推进 + y 出流外插 ... % ... x 方向推进 + x 出流外插 ... else % 偶数步：先 x 后 y % ... x 方向推进 + x 出流外插 ... % ... y 方向推进 + y 出流外插 ... end

这几十行看起来是重复劳动，但少了它，长时程积分的结果会慢慢歪掉。分裂法能不能用在长时间模拟里，很大程度上取决于这个交替有没有做到位。

5. 排查：五个把结果画飞的坑

5.1 图上只剩几根尖刺，其余是平板

现象：γ 取 1 或 2 时，surf 图上有几根极高的柱子，其余区域贴着零平面，看起来像数据缺了一块。

原因：放大因子对某些波数大于 1，这些模态指数增长，几步之内就把数值撑到极大。z 轴自动缩放后被最大值主导，其他结构全部压平，视觉上就成了平板。这不是缺数据，是尺度被吃了。

解决：把 axis 的 z 范围固定住，比如 axis([-20 20 -20 20 0 1])，让发散的部分冲出画面而不是压扁其余部分。同时在内层循环里打印 max(abs(U2(:)))，一旦这个值单调上升就是发散信号，比看图快得多。

5.2 波传到右上角又弹回来

现象：波形走到计算域边界后没有干净地离开，而是在角点附近堆积出一圈虚假的回波，随时间越来越明显。

原因：出流边界被当成了入流边界处理，给了固定值（通常是 0），等价于在出口立了一堵墙。或者虽然用了外插，但外插是在内点循环中途做的，用到了已经被覆盖的邻居值。

解决：出流边界一律用线性外插 2u_i - u_{i-1}，且必须放在内点循环全部结束之后，一次性覆盖整条边界。分裂法里两个方向的外插要分别在各自那一步之后做，不能合并到最后统一处理。

5.3 A 格式和 C 格式共用同一套循环上下界

现象：把 A 格式的公式替换进 C 格式的循环里，程序不报错，但边界附近出现异常值，或者结果从第一步就不对。

原因：C 格式只用到 i-1 和 j-1，内点范围是 2:end；A 格式和 Lax-Wendroff 用到 i+1 和 j+1，内点范围必须是 2:end-1，最后一圈留给边界处理。共用 2:end 会让 A 格式在 i = end 处访问越界元素，Matlab 会静默返回 0 或者报错，取决于写法。

解决：三种格式的循环范围分开写，中心类格式统一用 2:size(X,1)-1，边界另外补。

5.4 方台位置偏了或者形状是长方形的

现象：初值明明设的是正方形平台，画出来却偏在角落，或者变成宽窄不一的长条。

原因：meshgrid 的维度约定没对齐。meshgrid(x, y) 生成的是 X(i,j) = x(j)、Y(i,j) = y(i)，第一维对应 y，第二维对应 x。如果按「第一维是 x」的直觉去写 U1(13:21, 13:21)，两个方向的索引实际上是一致的所以看不出来；但只要 x、y 的网格点数不同，立刻错位。

解决：在脚本开头写一句 size(X) 和 X(1,:)、Y(:,1) 打印出来确认维度，比事后调图省事。统一约定「第一维 y、第二维 x」之后，所有索引都按这个来。

5.5 视频播出来抖动，或者根本只有一帧

现象：用 getframe 和 movie 做动画，窗口大小每一帧都在变，播出来画面跳；或者 writeVideo 写出的文件打不开。

原因：axis 范围随每帧数据自动缩放，figure 尺寸也跟着变，getframe 截到的帧大小不一致，播放时就会抖。视频文件则需要先 open(v) 再逐帧 writeVideo(v, m)，最后 close(v)，漏掉任何一步文件都是坏的。

解决：在循环外固定 set(gcf, 'Position', [...]) 和 axis([...])，让每一帧的尺寸完全一致；视频对象在循环开始前创建并 open，循环结束后 close，中间用 getframe 抓帧而不是抓图。

6. 取舍：精度、稳定域、实现成本，三者只能挑两个

四种格式放在一起，其实就是一张取舍表。迎风 C 格式稳定域最窄（γ ≤ 0.5），但只有一阶精度、耗散大、实现最简单；A 格式二阶精度、实现同样简单，但稳定域为零，根本不能用；Lax-Wendroff 二阶精度，稳定域进一步收紧到 0.3536，每步代价也翻倍；时间分裂法稳定域最宽（γ ≤ 1），每步可以用更大的步长，但时间精度掉回一阶，还多了一个方向交替的实现负担。

判断一个格式在具体问题上能不能用，光靠看表不够，得有一个能自动跑的监控量。稳定格式的总变差、L2 模、最大值都应该是有界的，波传播问题时总变差甚至会缓慢下降（耗散）；发散的格式这些量会单调指数上升。把下面这几行塞进时间循环，比事后看图靠谱得多：

tv = sum(sum(abs(diff(U2,1,2)))) + sum(sum(abs(diff(U2,1,1)))); l2 = sqrt(sum(U2(:).^2)); fprintf('n=%3d max=%.4e TV=%.4e L2=%.4e\n', n, max(U2(:)), tv, l2);

max 看有没有炸，TV 看波形有没有长出非物理的振荡，L2 看整体能量是不是守恒范围内的有界量。三个量一起看，稳定和不稳定几乎一眼就能分开，不用等到第 90 步画完图才发现白跑。

我自己的习惯是，每次换格式、换 γ、换初值，都先把这三行打印开着跑一遍，确认有界之后再关掉做动画。早年省这一步，用 γ = 1 跑 Lax-Wendroff，图上看着像「波传得挺整齐」，其实是数值已经长到 10 的十几次方、surf 的自动缩放把尖刺裁掉了——那次之后我就养成了先看数值再看图的毛病，改不掉了。希望这份拆解和代码能帮你少走一段同样的弯路，把源码里的循环和边界处理对着自己的推导跑一遍，比读十遍公式管用。

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