企业官网建设流程全解析

简介：这是一份面向计算流体力学初学者与高校课程复习者的Matlab上机实习报告，围绕二维对流方程在迎风C格式、A格式、Lax-Wendroff格式及时间分裂法下的数值求解展开，帮助读者理解差分格式构造、边界条件处理与稳定性判定。整份资料仅含1个PDF文件，压缩包约25.94MB，内含问题描述、初值条件、实验原理、运行结果、结果分析与附录程序等完整章节，可直接对照复现。作者按Δx=Δy=1、γ取0.5、1、2分别计算90个时间步长，给出四种格式的u-x-y波形图与动画素材，并借助误差传播方程推导放大因子，解释C格式仅在γ=0.5时稳定、其余发散的结论。附带的Matlab程序设计与实现细节，适合作为CFD课程作业参考、格式对比笔记与编程排错依据，已有319人学习。

1. 计算流体力学里的一个反直觉：同一个 γ，两种命运

计算流体力学（CFD）里最像玄学的一件事，是同一个参数 γ = Δt/Δx 在不同格式下判若两人。做二维对流方程 ∂u/∂t + ∂u/∂x + ∂u/∂y = 0 的数值解时，迎风 C 格式在 γ=1 已经炸成一片，时间分裂法却在 γ=1 给出几乎不衰减的方波平移；而 A 格式无论 γ 取 0.5 还是 2，理论上一开始就不该指望它收敛。同一道题、同一张 Matlab 网格，结论完全相反，原因不在代码写得好不好，而在格式本身的放大因子长什么样。

这张地图大致是这样：先用差分把偏导换成网格上的加减法，再用 von Neumann 放大因子判断这个换法靠不靠谱，最后用边界条件和可视化把结论坐实。它适合两类人：刚做完 CFD 上机、只会照着公式敲代码却说不清为什么发散的；以及复习计算流体力学、想把四种格式串成一条线的。

2. 方程与离散骨架：33×33 网格上的一块方波

二维对流方程的四种差分格式看起来公式一大堆，其实骨架只有三层：控制方程给物理，网格和初值给舞台，边界条件给“从哪进、从哪出”。这三层任何一层含糊，后面的稳定性分析就是在解释一个错误实验。梳理下来会发现，真正决定你能不能复现出报告里那几张曲面的，不是公式背得准不准，而是索引、步长换算和边界位置有没有对齐。

2.1 控制方程、初值与网格：把偏导搬到格点上

要算的方程是标量二维对流方程，两个方向的波速都取 1，也就是信息沿 x 正方向和 y 正方向同时传播：

∂u/∂t + ∂u/∂x + ∂u/∂y = 0，−16 ≤ x ≤ 16，−16 ≤ y ≤ 16，t > 0

初值是一个位于中心的方波：在 −4 ≤ x ≤ 4、−4 ≤ y ≤ 4 内 u=1，其余为 0。空间步长固定为 Δx = Δy = 1，于是每条边上有 33 个节点，网格是 33×33；方波占据的索引区间是 13:21，因为从 −16 数起，第 13 个节点正好落在 −4 上。这个换算关系一定要自己数一遍，索引错一位，后面所有图都是错的。

要特别注意 γ 的定义。报告里把 γ 记成 Δt/Δx = Δt/Δy，由于 Δx = Δy = 1，γ 在数值上就等于时间步长 Δt。所以“γ 取 0.5、1、2”实际是说时间步长取 0.5、1、2，考试里如果 Δx 不是 1，就绝不能把 γ 和时间步长混为一谈，这是最常见的换算翻车点。计算总长按要求跑到 90 个时间步，取帧和绘图频率则另算。

提示：Δx 固定为 1 只是教学网格的方便选择，稳定性判据里的 γ 永远是 Δt/Δx 这个无量纲比值，换网格后必须重新算。

2.2 从偏导到差分：把更新式看成一张账单

四种格式的差别，本质上是对 ∂u/∂x 和 ∂u/∂y 用几个点去近似。迎风 C 格式只取上游一点，A 格式取左右两点求中心差，Lax-Wendroff 在中心差基础上再加一个二阶修正项，时间分裂法则是把二维算子拆成两个一维算子分半步走。写成代码，循环骨架完全一样，只有中间那行更新式不同：

dx = 1; dy = 1; % 空间步长，网格共 33×33 个节点 r = 0.5; % r 即 gamma = dt/dx = dt/dy x = -16:dx:16; [X,Y] = meshgrid(x); % 生成二维坐标网格 U1 = zeros(size(X)); % n 时刻的解 U2 = zeros(size(X)); % n+1 时刻的解 U1(13:21,13:21) = 1; % 初值方波：x,y 均落在 [-4,4] for n = 1:90 for i = 2:size(U1,1) % 从 2 开始：第 1 行/列是入流边界 for j = 2:size(U1,2) U2(i,j) = U1(i,j) - r*(U1(i,j)-U1(i-1,j)) ... - r*(U1(i,j)-U1(i,j-1)); % 迎风 C 格式 end end U2(:,1) = 0; U2(1,:) = 0; % 入流边界：左列、下行给定为 0 U1 = U2; % 层间推进，U2 整体覆盖 U1 end

这段代码里三个参数值得逐个确认：size(U1,1)是行数（对应 y），size(U1,2)是列数（对应 x），U1(13:21,13:21)的先行后列顺序不能反；r一旦改成 1 或 2，同一段代码就跑出完全不同的结果。U1 = U2是整块替换而不是累加，写错成U1 = U1 + U2会得到一条指数增长的曲线，看起来像“格式发散”，其实是内存复用写错了。

2.3 边界条件：入流给值，出流外插

两个方向的波速都大于 0，意味着信息只从左边界和下边界进入计算域，从右边界和上边界离开。所以左边界（第 1 列）和下边界（第 1 行）是入流边界，直接给 u = 0；右边界和上边界是出流边界，不能硬给固定值，否则波打上去会反射回来，把内部的解污染掉。

迎风 C 格式只用到 i−1、j−1，出流边界天然不需要额外处理；但 A 格式和 Lax-Wendroff 格式用到了 i+1、j+1，必须在出流侧补一个条件，报告里用的是线性外插：

u(i+1,j) = 2u(i,j) − u(i−1,j)，u(i,j+1) = 2u(i,j) − u(i,j−1)

它的含义是“边界外面那一层，按内部最后两点的斜率线性延伸出去”，等于告诉格式：波继续往外走，不要回头。这段逻辑写成代码只有两行，但它是 A 格式和 Lax-Wendroff 能不能跑出“边界干净”的曲面的关键。

U2(:,end) = 2*U2(:,end-1) - U2(:,end-2); % 右边界（x 方向出流）线性外插 U2(end,:) = 2*U2(end-1,:) - U2(end-2,:); % 上边界（y 方向出流）线性外插 U2(:,1) = 0; % 左边界入流 U2(1,:) = 0; % 下边界入流

顺序上建议先算外插、再置入流，否则外插会把刚置零的入流点当成内部点参与计算。这两组语句谁先谁后，在 γ 较小时差别肉眼看不出来，等 γ 接近稳定上限时就会变成边界上一条明显翘起来的边。

3. 四种格式对比：迎风 C、A、Lax-Wendroff 与时间分裂

同一道题摆四种格式，比较的维度其实只有三个：精度阶数、数值耗散和稳定域。精度决定波峰会不会被削平，耗散决定削得多快，稳定域决定你能不能活着算到第 90 步。中间这两章的重点不在于把公式抄全，而在于看到一条公式就能预判它会画出什么形状的曲面——这件事在考场上比会推导更值钱，在机房调试时更是救命。

3.1 迎风 C 格式：用一阶精度换来的安全感

迎风 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)]

它做的事情非常朴素：把 x 方向的增量和 y 方向的增量分别扣掉。因为只取上游一点，它天然“顺着波走”，稳定条件是四个格式里最宽松的之一：γ ≤ 0.5。

代价是一阶精度带来的强数值耗散。方波的四个角会被磨圆，波峰高度随时间下降，看起来像被热传导抹过一遍。用修改方程的话说，一维迎风等价于在原方程右边加了一项 (1−γ)aΔx/2 · u_xx 的数值粘性，γ 越接近 1，这个粘性越小；γ = 1 时一维迎风是精确平移，一点不糊。但二维迎风 C 是两个方向耗散的叠加，稳定上限被压到 0.5，所以二维里根本到不了“零耗散”那个甜点，这也是为什么报告里 γ=0.5 的曲面已经明显掉高度、边界也不再锐利。

复现时最该盯住的现象是：γ=0.5 时方波整体向右上平移、峰高缓慢下降、四个角变钝；γ=1 和 γ=2 时数值迅速抬高并出现正负振荡，几十步之内就冲出坐标范围。

3.2 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)]

从截断误差看它是二阶精度，比迎风 C 漂亮得多，对称、无偏、不引入方向性耗散。但实验里 γ 取 0.5、1、2 三个值全部发散，这不是调参能救的，而是格式本身的缺陷。

原因在放大因子上：G = 1 − iγ(sinβ + sinθ)，取模得 |G| = sqrt(1 + γ²(sinβ+sinθ)²)，只要 γ>0 且波数不为零，|G| 就严格大于 1。也就是说中心差分对任何非零波数都在放大误差，网格越细、初值越光滑，它照样放大，只是放大得慢一点。这类格式在数值上毫无“后悔药”可吃。

一个很实用的观察点：A 格式跑出来的图不是简单地“变成 NaN”，而是先出现高频锯齿状振荡，再整体抬高。如果只看前几步的图觉得“还行”，很容易误判成稳定，等跑到第 40 步才发现已经飞了。

3.3 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 格式稳。但二维的稳定域比一维窄得多。按 K1 = Δt/Δx、K2 = Δt/Δy 的定义，二维 Lax-Wendroff 的判据是：

K1^(2/3) + K2^(2/3) ≤ 1

当 K1 = K2 = γ 时，条件变成 2γ^(2/3) ≤ 1，解出 γ ≤ 2^(−3/2) ≈ 0.3536。这就解释了为什么 γ=0.5、1、2 三个值全部发散——它们都在这条线之上，而补做的 γ=0.3 实验就收敛了。

注意：一维 Lax-Wendroff 的稳定上限是 γ ≤ 1，把同一套代码直接搬到二维，稳定域会缩到 0.3536，这是最容易踩的坑：代码没改错，是判据换了。

收敛之后也不能高兴太早。Lax-Wendroff 没有迎风那种耗散，方波边缘会出现明显的色散波纹，峰前后各拖出一串小幅振荡，这是二阶格式的典型痕迹，不是程序写错了。

3.4 时间分裂法：二维拆成两个一维，γ=1 反而最准

时间分裂法把二维问题拆成两步一维问题，每个半步只处理一个方向：

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)]

代码上它多了一个中间变量 U3 存放半时间步的解，每步先算一个方向，用半时间步的结果再算另一个方向。实验里 γ=0.5 和 γ=1 都收敛，γ=2 发散，稳定条件是 γ ≤ 1。

最反直觉的地方在这里：γ=1 时，一维迎风格式恰好退化成精确平移 u(n+1,i) = u(n,i−1)，耗散为零；两个方向各平移一格，整体就是把方波沿对角线挪一步，形状几乎不变。所以时间分裂法在 γ=1 反而给出了四种格式里最干净的曲面，而不是像迎风 C 那样在 γ=1 直接炸掉。这条结论几乎每年都会有人在复习时记反。

还有一个实现细节值得抄进笔记：为了避免固定按“先 y 后 x”的顺序造成方向上的系统偏差，程序在奇数步用 y→x 顺序、偶数步用 x→y 顺序。这一步不影响稳定性判据（放大因子是两个一维因子相乘，交换顺序不改变模），但会影响图形的对称性，尤其是方波的对角线走向。

4. 稳定性分析：von Neumann 放大因子怎么一步步算

稳定性这件事，实验图像只能告诉你“炸了”或“没炸”，说不清临界值在哪。von Neumann 方法的价值在于：它把“会不会炸”变成一个可以手算的不等式，让你在写代码之前就知道 γ 该取多少，而不是跑完 90 步再回头猜。的推导套路只有三步，四种格式走的是同一条流水线，区别只在代入后整理出的 G 长什么样。

4.1 误差传播方程与三步法

设数值解和真解之间的误差为 ε，由于方程是线性的，误差本身满足和原差分格式完全相同的方程。以迎风 C 格式为例，误差传播式就是把 u 换成 ε 的那条式子：

ε(n+1,i,j) = ε(n,i,j) − γ[ε(n,i,j) − ε(n,i−1,j)] − γ[ε(n,i,j) − ε(n,i,j−1)]

第一步，设误差按单一波数传播：ε(n,i,j) = A^n · exp[i(k1·x + k2·y)]。第二步，把相邻点的误差用同一个 A^n 表示，例如 ε(n,i−1,j) = A^n·exp[i(k1(x−Δx) + k2·y)]，这样所有项都能约掉空间因子。第三步，剩下的是层间放大因子 G = A^(n+1)/A^n，它是一个复数，稳定性判据就是 |G| ≤ 1 对所有波数成立。

记 β = k1Δx、θ = k2Δy，把指数用欧拉公式拆成实部和虚部，迎风 C 的放大因子整理成：

G = 1 − 2γ + γ(e^(−iβ) + e^(−iθ))

这一步是整个推导里最容易写错的地方：符号漏一个负号，后面所有结论都会反。建议整理完先代两个极端波数自检——β = θ = 0（零波数，应该 |G| = 1）和 β = θ = π（最细锯齿波，通常给出最坏情况）。

想验证自己的推导对不对，不用肉眼看公式，直接拿几行 Matlab 把 |G| 的最大值扫出来就行：

beta = linspace(-pi, pi, 201); theta = linspace(-pi, pi, 201); [B,T] = meshgrid(beta, theta); % 覆盖所有波数组合 glist = 0:0.05:1.2; for k = 1:numel(glist) r = glist(k); G = 1 - r*(1-exp(-1i*B)) - r*(1-exp(-1i*T)); % 迎风 C 的放大因子 amax(k) = max(abs(G(:))); % 全波数上的最大模 end plot(glist, amax, glist, ones(size(glist))); % 与 |G|=1 的水平线比较

图上曲线第一次触到 1 的那个 γ，就是稳定上限，迎风 C 会给出 0.5。把G = ...那一行换成 A 格式的1 - 1i*r*(sin(B)+sin(T))，曲线会整体在 1 之上，直观看到“无条件不稳定”。

4.2 四个格式的 |G| 与稳定阈值

把四种格式的放大因子并排放，稳定性的差别一眼就出来了：

迎风 C 的临界值这样得到：取最坏的锯齿波 β = θ = π，G = 1 − 4γ，要求 |1 − 4γ| ≤ 1，得 0 ≤ γ ≤ 0.5。γ = 0.5 时进一步可验证 |G|² = 0.5[1 + cos(β−θ)] ≤ 1 对所有波数成立，所以 0.5 是既不充分条件里的紧界，理论值与实验完全对上。

A 格式的结论更干脆：|G|² = 1 + γ²(sinβ+sinθ)²，只要 γ>0，总能取到 sinβ+sinθ ≠ 0，模严格大于 1，无条件不稳定。这也解释了为什么三个 γ 值全都发散——不是取值不够小，是根本没有安全区。

时间分裂法要分开看两个因子：一维迎风因子的模平方是 1 − 2γ(1−γ)(1−cosβ)，只要 0 ≤ γ ≤ 1 就不超过 1；两个因子相乘，模也不会超过 1。所以 γ ≤ 1 是它的稳定界，γ=2 时单个因子的模已经大于 1，乘起来发散，与实验一致。顺带解释了 γ=1 的“零耗散”：此时 1−cosβ 的系数 2γ(1−γ) 恰好为零，|G| 恒等于 1。

4.3 γ=0.5/1/2 的实验现象与理论对照

把报告里的图像现象和理论判据并排放，是复习时最省力的一张表：

这张表里有两处最容易记反：Lax-Wendroff 在 γ=0.5 是发散的，虽然 0.5 听上去“不大”；时间分裂法在 γ=1 是稳定的，而且是最准的一档。判断依据永远回到 |G| ≤ 1，而不是“γ 越小越安全”这种直觉——γ 太小反而会放大迎风格式的耗散，让方波糊得更快。

5. 易漏点：实验报告里最容易翻车的五个地方

同一份代码，有人跑出干净的曲面，有人跑出一团噪声，差别往往不在格式，而在这五个细节。它们有个共同特点：程序不报错、Matlab 不警告、图还能画出来，所以极难自查。下面按“错在哪 → 为什么容易漏 → 怎么补”逐条过，每条都是可以直接照着改的。

5.1 入流边界到底该置哪一行

易漏点：把U2(end,:) = 0当成入流边界条件写上，而end对应的是 y 最大的那一行，也就是出流侧。

为什么漏：入流和出流在公式里只是“左/右”“上/下”的区别，写代码时1和end长得同样顺手；而且置零之后程序照跑，图也照出，只是边界附近多了一条奇怪的平台，容易被当成数值耗散。

怎么补：波速为正时，入流永远是第 1 列和第 1 行，写成U2(:,1)=0; U2(1,:)=0;，出流永远是end那一侧，用线性外插。改完后对着 γ=0.5 的曲面看四个边，入流侧应该是一条干干净净的零线，出流侧应该看到波“走出去”而不是堆在边上。

5.2 出流外插的 end-2 写成 end-1

易漏点：外插式写成2*U2(end-1,:) - U2(end-1,:)，两个点写成同一个点。

为什么漏：这样写等价于U2(end-1,:)，也就是“边界值等于里面一格”，零梯度边界，看上去也合理，甚至还更稳；只有在波传到边界时才会发现反射比预期明显，但那时候已经很难判断是格式问题还是边界问题。

怎么补：严格按 u(i+1) = 2u(i) − u(i−1) 写，第二项必须是end-2，两个方向的边界各检查一遍。一个笨但有效的自查方法：临时把外插换成U2(end,:) = U2(end-1,:)跑一次，如果两张图几乎一样，说明你的波还没走到边界，这时候不能算“边界条件验证通过”。

5.3 循环从 2 开始，却漏掉第 1 行的更新

易漏点：内层循环写for i = 2:size(U1,1)，于是 i=1 和 j=1 的整行整列永远保持初值或上一层的旧值，忘了每步都要重新赋 0。

为什么漏：初值里这两行本来就是 0，前几十步看不出异常；但入流边界旁边那一圈如果不每步刷新，等于边界条件只在第 0 步生效，后面全靠旧数据撑着。

怎么补：把边界赋值放在时间循环内部、内层循环之后，每步都执行一次，顺序上先外插后置零。想确认有没有漏，把U2(:,1)单独打印出来看是否恒为 0，比盯着曲面图可靠得多。

5.4 时间分裂法的半时间步变量与奇偶顺序

易漏点：两个半步共用一个数组，或者忘了在奇数步和偶数步交换 x/y 的先后顺序。

为什么漏：半步的结果如果直接覆盖原数组，第二步就会拿已经被修改过的邻居去算，等于把两个方向耦合在一起，格式悄悄变成了另一种东西，稳定性判据也就不再适用。

怎么补：用三个数组分工——U1 是整步初值，U3 存半时间步结果，U2 存整步结果，每步结束再U1 = U2。顺序上用mod(n,2)判断奇偶，奇数步先 y 后 x、偶数步先 x 后 y，并在两个半步各自补上出流外插。这样跑出来的 γ=1 结果应该是一个几乎不失真的对角平移方波，如果峰高掉得厉害，先查这一步。

5.5 把 0.3536 当成拍脑袋的数

易漏点：Lax-Wendroff 直接用一维的 γ ≤ 1 去设参数，跑到 γ=0.5 发现发散，怀疑代码有 bug。

为什么漏：一维和二维的稳定域不是同一个判据，二维多了一个方向的波数耦合，条件从线性收紧成 2/3 次幂形式，数值上从 1 掉到 0.3536，差了三倍。

怎么补：把判据 K1^(2/3) + K2^(2/3) ≤ 1 写在代码注释里，等步长时直接记 γ ≤ 2^(−3/2) ≈ 0.3536；调参时先取 0.3 确认格式本身没问题，再往临界值靠。凡是“理论说稳、实验说炸”的情况，第一反应都应该是回去核对判据的维度，而不是重写循环。

6. 复习清单：从上机报告到考场推导的自查路径

把这份材料当复习提纲用，最有效的顺序不是从头读公式，而是拿一张白纸，按“方程 → 网格 → 格式 → 放大因子 → 临界值 → 现象”六步默写一遍，每一步都问自己能不能说出下一层为什么这么定。下面这张表是我自己复盘时用的检查项，左边是动作，中间是判据，右边是出错时的典型症状，对着走一遍基本能定位到问题在哪一层。

默写完公式，再用代码验一遍，比纯看公式记得牢。最小验证脚本不用长，把放大因子扫一遍、把边界打印出来看，两件事就够：

% 1) 扫出迎风 C 的稳定上限，应停在 0.5 附近 r = 0.5; [B,T] = meshgrid(linspace(-pi,pi,101)); G = 1 - r*(1-exp(-1i*B)) - r*(1-exp(-1i*T)); fprintf('gamma=%.2f, max|G|=%.4f\n', r, max(abs(G(:)))); % 2) 确认入流边界每步都被刷新 assert(all(U2(:,1) == 0) && all(U2(1,:) == 0), '入流边界没有每步置零');

第一段把 γ 换成 0.6、0.8、1.0 各跑一次，max|G| 会依次越过 1，临界点就“长”在手上了；第二段的断言放在时间循环里，能第一时间抓住 5.3 里那个“只在第 0 步生效”的边界。这两件事加起来不超过二十行，却比反复调参省掉一晚上。

还有一个查漏补缺的笨办法：把四种格式的曲面按 γ=0.5、1、2 摆成 3×4 的九宫格（Lax-Wendroff 补一张 γ=0.3），不看标题只看形状，先自己判断哪个稳、哪个炸，再翻回第 4 章对答案。能全部判对，说明放大因子那套推导是真的懂了，而不是背下了四行结论。

我自己的教训是：第一次做这个上机时，把出流外插的end-2写成了end-1，图看着挺正常，直到答辩前被问“边界为什么有反射”才翻出来，白跑了两天的动画。从那以后我每次写二维循环，都强制先把四个边界各打印一行、确认每步都在刷新，再开始调 γ。希望帮到你。

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