One Forth

P11 · 第 1 章 · 第一编 数值建模

离散化、误差预算与模型验证

把物理模型到数值结论的误差链逐层拆开，用一致性、稳定性、网格加密、守恒残差、制造解与独立基准验证实现，并区分数值 verification 与面向实验的 validation。

难度 5/5约 220 分钟最近修订 2026-07-14
报告页面错误
预备知识稳定性、收敛性与科学计算综合复习理想化、守恒量与模型边界浮点数、条件数与误差传播

本章目标

  1. 区分建模、参数、离散、迭代、舍入和统计误差，并把实现缺陷与可接受误差分开。
  2. 用 Taylor 展开解释一致性，以扰动放大解释稳定性，并在适用前提下连接到收敛。
  3. 设计至少三层网格或时间步加密，计算观测阶并使用 Richardson 外推核对渐近区。
  4. 用量纲一致的守恒残差、制造解和解析基准验证代码与单次计算结果。
  5. 区分 code/solution verification 与 validation，建立带单位、相关性和适用范围的误差预算。
页面阅读位置0% · 仅保存在此浏览器
章节未开始
本册完成进度0/6 章 · 0%
本页目录

一条数值结论有多层来源

计算物理不是把公式输入程序后得到一个“精确答案”。实验对象先被理想化为物理模型，再写成连续或随机方程，随后经过空间与时间离散、代数或非线性迭代、有限精度运算及统计汇总。每层回答不同问题：

  1. 建模误差来自忽略的物理、闭合关系和适用范围。
  2. 参数与初边值不确定度来自测量、标定或估计。
  3. 截断或离散误差来自用有限网格、时间步和有限基逼近连续问题。
  4. 迭代误差来自在线性、非线性或优化求解器未完全收敛时停止。
  5. 舍入误差来自有限浮点表示、消去和归约次序。
  6. 统计误差来自有限随机样本或有限相关轨迹。

实现错误不是应当写进预算后接受的第七类误差。数组越界、边界符号写反或单位转换错误应通过测试和 verification 找出并修复。误差预算描述在已知实现下仍不可避免或有意保留的近似。

若目标量 QQ 的单位是 W m−2\mathrm{W\,m^{-2}}，绝对误差 Qh−QQ_h-Q 也必须是 W m−2\mathrm{W\,m^{-2}}；相对误差 ∣Qh−Q∣/Qref|Q_h-Q|/Q_{\mathrm{ref}} 才无量纲。方程残差常带自己的单位，直接比较不同变量的原始残差没有意义，应按特征尺度或容差归一。

Taylor 截断误差与一致性

以光滑函数的一阶导数为例，中心差分

Dhf(x)=f(x+h)−f(x−h)2hD_hf(x)=\frac{f(x+h)-f(x-h)}{2h}

由 Taylor 展开得到

Dhf=f′(x)+h26f(3)(x)+O(h4).D_hf=f'(x)+\frac{h^2}{6}f^{(3)}(x)+O(h^4).

若 x,hx,h 以米计而 ff 以开尔文计，DhfD_hf 和截断误差单位均为 K m−1\mathrm{K\,m^{-1}}。称格式对微分算子一致，是指把足够光滑的精确解代入离散方程后，局部截断误差随 h,Δt→0h,\Delta t\to0 消失。一致性说明方程逼近对了对象，却不单独保证多步计算中的扰动不被放大。

例 1：用两次加密观测中心差分阶数

取无量纲变量 f(x)=sin⁡xf(x)=\sin x，在 x=0.5x=0.5 求导。精确值 f′(0.5)=cos⁡0.5≈0.877583f'(0.5)=\cos0.5\approx0.877583。当 h=0.10h=0.10，

D0.10f=sin⁡0.60−sin⁡0.400.20≈0.876121,D_{0.10}f =\frac{\sin0.60-\sin0.40}{0.20} \approx0.876121,

绝对误差约 1.46×10−31.46\times10^{-3}。当 h=0.05h=0.05，

D0.05f=sin⁡0.55−sin⁡0.450.10≈0.877217,D_{0.05}f =\frac{\sin0.55-\sin0.45}{0.10} \approx0.877217,

误差约 3.66×10−43.66\times10^{-4}。误差比约为 3.993.99，故观测阶 pobs=log⁡2(3.99)≈2.00p_{\mathrm{obs}}=\log_2(3.99)\approx2.00，与二阶截断项一致。若继续减小 hh 后误差不再按四倍下降，可能已离开截断误差主导区，不能仍宣称观测到二阶。

稳定性、收敛与适用前提

稳定性研究离散推进对初值误差、边界误差和每步舍入扰动的放大。对线性、适定的初值问题，Lax 等价框架说明一致且稳定的线性差分格式收敛，反之亦然。这个结论有明确前提，不能直接替代非线性方程、移动间断、随机算法或非适定逆问题的专门分析。

考虑一维热方程

∂u∂t=α∂2u∂x2,[α]=m2 s−1.\frac{\partial u}{\partial t} =\alpha\frac{\partial^2u}{\partial x^2}, \qquad [\alpha]=\mathrm{m^2\,s^{-1}}.

显式中心格式

uin+1=uin+r(ui+1n−2uin+ui−1n),r=αΔtΔx2u_i^{n+1} =u_i^n+r(u_{i+1}^n-2u_i^n+u_{i-1}^n), \qquad r=\frac{\alpha\Delta t}{\Delta x^2}

在一维周期或相容边界下要求 0≤r≤1/20\le r\le1/2 才满足 Fourier 稳定条件。rr 无量纲。稳定上限只是防止某些误差模指数放大，不代表接近上限时目标量已足够准确。

例 2：热扩散显式步长上限

取 α=1.00×10−5 m2 s−1\alpha=1.00\times10^{-5}\,\mathrm{m^2\,s^{-1}}、 Δx=1.00 mm=10−3 m\Delta x=1.00\,\mathrm{mm}=10^{-3}\,\mathrm m。条件给

Δt≤Δx22α=10−62.00×10−5=0.0500 s.\Delta t\le \frac{\Delta x^2}{2\alpha} =\frac{10^{-6}}{2.00\times10^{-5}} =0.0500\,\mathrm s.

选 Δt=0.0400 s\Delta t=0.0400\,\mathrm s 时 r=0.400r=0.400，满足线性稳定性。若网格加密为 Δx/2\Delta x/2，保持同一 rr 必须把时间步减为原来的四分之一；只加密空间而不调整时间步会破坏稳定条件，也无法把空间误差与时间误差分开。

三层加密与观测收敛阶

设某标量输出在渐近区满足

Qh=Q+Chp+O(hp+1).Q_h=Q+C h^p+O(h^{p+1}).

对比例为 2 的三层网格 h,h/2,h/4h,h/2,h/4，

pobs=ln⁡∣(Qh−Qh/2)/(Qh/2−Qh/4)∣ln⁡2.p_{\mathrm{obs}} =\frac{\ln\left| (Q_h-Q_{h/2})/(Q_{h/2}-Q_{h/4}) \right|}{\ln2}.

至少三层才能同时观察阶数与细化趋势。若差值变号、阶数剧烈变化或局部自适应网格不能用单一 hh 表示，应报告更完整的误差范数与自由度，而不是强行套一个整数阶。比较网格时还要保持物理参数、边界、求解容差和后处理定义一致。

场误差还必须声明范数。对控制体误差 eie_i 和体积 ViV_i，一种体积加权均方范数是

∥e∥2,V=(∑iei2Vi∑iVi)1/2,∥e∥∞=max⁡i∣ei∣.\|e\|_{2,V} =\left( \frac{\sum_i e_i^2V_i}{\sum_iV_i} \right)^{1/2}, \qquad \|e\|_\infty=\max_i|e_i|.

二者单位都与被比较场相同。直接对网格向量用未加权 Euclidean 范数，会随节点数改变尺度，跨网格不一定可比。L2L^2 范数强调总体误差，L∞L^\infty 对局部峰值敏感；含间断的解在不同范数中可能呈现不同阶数。

一个积分输出偶然达到高阶，也不能证明整个场同阶收敛。反之，局部尖角附近点值阶数降低，某个守恒通量仍可能准确。验证计划应同时选与决策相关的目标量、至少一个全场范数和守恒残差，并预先写明接受阈值。

例 3：热流的观测阶与 Richardson 外推

某壁面平均热流在三层网格上为

Qh=1.200,Qh/2=1.050,Qh/4=1.0125kW m−2.Q_h=1.200,\quad Q_{h/2}=1.050,\quad Q_{h/4}=1.0125 \quad \mathrm{kW\,m^{-2}}.

相邻差值为 0.1500.150 与 0.0375 kW m−20.0375\,\mathrm{kW\,m^{-2}}，比值 4，故 pobs=2p_{\mathrm{obs}}=2。用最细两层作外推：

Qext=Qh/4+Qh/4−Qh/22p−1=1.0000 kW m−2.Q_{\mathrm{ext}} =Q_{h/4} +\frac{Q_{h/4}-Q_{h/2}}{2^p-1} =1.0000\,\mathrm{kW\,m^{-2}}.

最细网格相对外推值高 1.25%1.25\%。这个估计依赖单一二阶项主导；三点恰好成二阶序列是证据，不是对所有更细网格的证明。

迭代误差、条件数与舍入下限

离散后常得到 Ax=bA\boldsymbol x=\boldsymbol b。残差 r=b−Ax^\boldsymbol r=\boldsymbol b-A\hat{\boldsymbol x} 小，不必保证解误差小；粗略关系受条件数 κ(A)\kappa(A) 放大。停止准则应使用与变量尺度相容的相对残差、物理残差或预条件范数，并通过更严容差复算确认目标量不再变化。若离散误差约为 10−310^{-3}，把迭代误差压到 10−1210^{-12} 可能浪费计算；反之，容差为 10−210^{-2} 会掩盖网格阶。

减小 hh 也不是无限改善。差分中的两个相近浮点数相减会放大相对舍入误差，导数误差常呈先下降后上升的 U 形。双精度、缩放、求和顺序和矩阵条件数都应记录。用更高精度或改变公式能诊断舍入主导，但不能修复错误的模型或边界条件。

守恒残差是带单位的独立检查

对有限体积或粒子系统，全局质量、动量、能量与电荷平衡提供不依赖局部点值的检查。固定控制体质量缺陷可写为

δm=M(t2)−M(t1)−∫t1t2(m˙in−m˙out) dt,\delta_m =M(t_2)-M(t_1) -\int_{t_1}^{t_2} (\dot m_{\mathrm{in}}-\dot m_{\mathrm{out}})\,\mathrm dt,

单位为千克。若按时间报告平衡残差 δm/(t2−t1)\delta_m/(t_2-t_1)，单位变为 kg s−1\mathrm{kg\,s^{-1}}。归一分母要声明，可选初始存量、累计通量或特征质量；分母接近零时应保留绝对残差。

例 4：储罐计算的质量平衡缺陷

初始质量 10.0 kg10.0\,\mathrm{kg}，四秒内恒定流入 2.00 kg s−12.00\,\mathrm{kg\,s^{-1}}、流出 0.500 kg s−10.500\,\mathrm{kg\,s^{-1}}。精确账本要求末质量

M⋆=10.0+(2.00−0.500)×4.00=16.0 kg.M_\star=10.0+(2.00-0.500)\times4.00 =16.0\,\mathrm{kg}.

某计算给 15.94 kg15.94\,\mathrm{kg}，则 δm=−0.060 kg\delta_m=-0.060\,\mathrm{kg}，平均缺陷率 −0.015 kg s−1-0.015\,\mathrm{kg\,s^{-1}}，相对末质量为 −0.375%-0.375\%。这项检查能发现净丢失，却不能仅凭缺陷位置断定是通量、时间推进还是边界实现错误，还需局部残差和加密实验定位。

守恒到机器精度也不等于解正确：一个符号写错但内部成对抵消的程序仍可能守恒。守恒、收敛阶和独立基准必须联合使用。

制造解、解析基准与 verification

制造解方法先任取足够光滑且满足所需边界类型的函数 um(x,t)u_m(\boldsymbol x,t)，再把它代入目标微分算子，反算源项

fm=L(um).f_m=\mathcal L(u_m).

数值代码用 fmf_m 和相应初边值求解，检查误差是否按设计阶下降。制造解可以覆盖复杂源项、非均匀系数和边界分支，不要求 umu_m 是原始无源物理问题的真实解。它验证“离散方程和代码是否按预期求解”，不验证模型能否描述自然。

解析解、独立高精度求解器、对称极限和已知守恒解也是基准。Code verification 关注实现与离散阶；solution verification 估计某次生产计算的网格、时间步和迭代不确定度。只在一个网格与一个参考数值吻合，无法区分偶然误差抵消。

Validation 的证据边界

Validation 比较模型预测与实验或观察数据，问的是“模型在指定用途和条件下是否足够”。它必须记录测量不确定度、参数校准数据与验证数据是否独立、比较量和单位、空间时间对齐以及接受标准。用同一数据既拟合参数又宣称独立验证，会高估证据。

verification 通过不能保证 validation 通过：方程可以被准确求解但遗漏关键物理。反过来，单个实验点吻合也不能证明程序正确，离散误差和模型偏差可能抵消。Validation 支持的是限定工况、输出和误差容限内的使用主张，不是对模型“永远正确”的认证。

误差预算与报告

每个误差项至少记录来源、估计方法、单位、符号或区间、是否随机以及相关性。只有当各项可视为独立零均值标准不确定度且单位一致时，才可用平方和开根号：

uc=u12+⋯+un2.u_c=\sqrt{u_1^2+\cdots+u_n^2}.

已知偏差应校正或单列，有界模型差应作为情景区间，不能都伪装成 Gaussian 标准差。

例 5：温度预测的分层误差预算

某温度输出为 T=350.0 KT=350.0\,\mathrm K。网格加密估计离散标准不确定度 0.40 K0.40\,\mathrm K，更严容差给迭代影响 0.05 K0.05\,\mathrm K，重复随机输入给统计标准误 0.30 K0.30\,\mathrm K，舍入诊断小于 0.01 K0.01\,\mathrm K。若前三项近独立且零均值，

unum≈0.402+0.052+0.302+0.012=0.502 K.u_{\mathrm{num}} \approx\sqrt{0.40^2+0.05^2+0.30^2+0.01^2} =0.502\,\mathrm K.

另有模型简化引起的有向偏差区间约 [−1.0,0.2] K[-1.0,0.2]\,\mathrm K，应单列而不并入上述平方和。报告可写成“数值标准不确定度约 0.50 K0.50\,\mathrm K，另有模型偏差情景区间”，比给一个来源不明的 ±1.1 K\pm1.1\,\mathrm K 更可审计。

常见误区

常见误区

“残差降到机器精度就表示物理解精确。”残差只说明离散代数方程被解得很紧，模型和离散误差仍可更大。

常见误区

“两个网格给相近结果已经证明收敛阶。”至少需要三层观察阶数，并确认处于渐近区且其他误差更小。

常见误区

“实验吻合就是 verification。”实验比较属于 validation；代码与离散实现需要制造解、解析基准和加密研究独立验证。

练习：从截断误差到验证主张

练习
所属知识
误差分类
难度
2/5

分别给“忽略辐射、有限网格、求解器提前停止、相近浮点数相减、有限样本均值”分类。

查看提示
判断误差发生在物理方程、网格、代数求解、浮点或采样哪一层。
查看解答
忽略辐射是建模误差；有限 Δx\Delta x 是离散误差；线性迭代提前停止是迭代误差；相消损失是舍入误差；有限随机样本均值波动是统计误差。
练习
所属知识
观测阶
难度
3/5

三层输出 1.120,1.040,1.020 Pa1.120,1.040,1.020\,\mathrm{Pa}，求观测阶和二阶 Richardson 外推值。

查看提示
用 p=ln⁡(∣Qh−Qh/2∣/∣Qh/2−Qh/4∣)/ln⁡2p=\ln(|Q_h-Q_h/2|/|Q_h/2-Q_h/4|)/\ln 2。
查看解答
差值为 0.080 和 0.020，比例 4，所以 p=2。外推值为 1.020+(1.020−1.040)/3=1.013331.020+(1.020-1.040)/3=1.01333，单位与 Q 相同。
练习
所属知识
扩散稳定性
难度
3/5

扩散率 2.0×10−6 m2 s−12.0\times10^{-6}\,\mathrm{m^2\,s^{-1}}、网格 0.50 mm0.50\,\mathrm{mm}，求显式时间步上限及网格减半后的上限。

查看提示
一维显式中心格式要求 αΔt/Δx2≤1/2\alpha \Delta t/\Delta x^{2}\le 1/2。
查看解答
α=2×10−6m2/s\alpha=2\times 10^{-6} m^{2}/s、Δx=0.5mm\Delta x=0.5 mm 时 Δt≤(2.5×10−7)/(4×10−6)=0.0625 s\Delta t\le(2.5\times 10^{-7})/(4\times 10^{-6})=0.0625\,\mathrm{s}。空间加密一半时上限变为 0.015625 s。
练习
所属知识
守恒单位
难度
3/5

写出质量平衡缺陷及缺陷率的单位，并说明净通量接近零时如何报告。

查看提示
存量变化和时间积分通量都是 kg；除以时间后才是 kg/s。
查看解答
δm=M2−M1−∫(m˙in−m˙out)dt\delta_m=M_2-M_1-\int(\dot{m}_{in}-\dot{m}_{out})dt，单位 kg；δm/Δt\delta_m/\Delta t 单位 kg/s。若累计净通量近零，不宜用它归一，应报告绝对缺陷或选非零特征质量。
练习
所属知识
制造解
难度
4/5

说明制造解的构造步骤，以及它为何不能替代 validation。

查看提示
先选 umu_m，再由目标算子计算源项与边界值；不要声称它验证自然模型。
查看解答
制造解检查算子、源项、边界和预期收敛阶，即 code verification；因为源项是为所选函数人工构造，它不证明原物理模型与实验一致。
练习
所属知识
误差预算
难度
4/5

给出把离散、迭代、测量和模型误差合并前必须检查的条件，并指出哪些项不宜直接平方和。

查看提示
先统一单位并区分独立标准不确定度、有向偏差和有界区间。
查看解答
独立零均值标准项可平方和开根号；已知偏差应校正或单列；模型区间和强相关项不能当作独立 Gaussian 项。报告需保留来源、单位、相关性和适用条件。

知识连接与资源

课程 · 2015

Numerical Fluid Mechanics

Pierre Lermusiaux

用于核对 P11 的离散化误差、场方程网格算法、时间步稳定性、边界条件和数值验证流程。

打开官方来源

MIT OpenCourseWare 2.29 覆盖数值误差、稳定性、网格、时间推进与验证，可用于核对本章离散分析和加密流程。本章示例使用明示教学参数；未给实验来源的数字不构成真实系统 validation 数据。