› Richardson 外推与观测收敛阶模拟器 返回
数值分析

Richardson 外推与观测收敛阶模拟器

从3个步长得到的数值解实时计算观测收敛阶 p 和 Richardson 外推真值推定 u_ext。用于 CFD、FEM 的网格收敛（V&V）、ODE 解法的精度评估，以及中点法则、梯形法则和 RK 法的阶数检验。

参数设置
粗步长 h₃
最粗的步长（h₃ = 最大）
中步长 h₂
中间步长（通常 h₃/2）
细步长 h₁
最细的步长（h₁ = 最小）
粗解 u₃
h₃ 得到的数值解
中解 u₂
h₂ 得到的数值解
细解 u₁
h₁ 得到的数值解
计算结果
—
步长比 r = h₂/h₁
—
观测收敛阶 p
—
Richardson 外推值 u_ext
—
细解截断误差
—
渐近误差比 e₁/e₂
—
阶数判定
—
log-log 图 — 误差 vs 步长 h（斜率 = p）

蓝点：3个步长 h₁,h₂,h₃ 的误差。黄线：通过3点的回归直线（斜率 = 观测收敛阶 p）。绿菱：Richardson 外推值（h→0）。

解 vs 步长 h（含外推值）
误差 vs 步长 h（双对数坐标）
理论与主要公式

$$u_{exact} \approx u_1 + \frac{u_1-u_2}{r^p-1},\qquad p=\frac{\ln\bigl(|u_2-u_3|/|u_1-u_2|\bigr)}{\ln r}$$

r 是步长比（h₂/h₁），p 是观测收敛阶。u₁,u₂,u₃ 分别为细、中、粗数值解。当3解处于渐近区域时，u_ext 用误差小 1~2 个数量级的精度近似真值。

$$e_h \approx C\,h^{p},\qquad \frac{e_1}{e_2}=\frac{1}{r^{p}}$$

离散误差 e_h 按 h 的 p 次方减小。渐近误差比 e₁/e₂ 接近 1/r^p 表示进入渐近区域，观测阶值可信。

$$\mathrm{GCI}=\frac{F_s\,|u_1-u_2|}{(r^p-1)\,|u_1|}$$

网格收敛指数（Roache 1998）。安全系数 F_s=1.25（3网格）给出 Richardson 外推误差的信任区间。CFD V&V 报告的标准。

Richardson 外推与收敛阶

🙋
"Richardson 外推"是什么意思？为什么要用3个不同的步长？1个不就够了吗？
🎓
好问题。一个步长 h 计算的数值解 u_h 包含未知的误差 e_h。要了解"有多大偏离"和"偏离方向"，至少需要2个步长来看解的差异。而要"把收敛阶 p 作为未知数来估算"就需要3个点。比如 h=0.1 得 u=0.8125，只从这一点无法判断真值是 0.81、0.82 还是 0.85。但如果有3个数据点，就能反推出"啊，这是以 p=2 的速度收敛到 0.811"。
🙋
我看到"观测收敛阶"和"理论收敛阶"两个概念。教科书说"中点法则是2阶精度"，那是理论阶吗？
🎓
完全正确。"理论阶"是纸上推导的值（中点法则=2、RK4=4 等），"观测阶"是实际跑代码的结果。两者一致的前提是：(1)程序无 bug，(2)步长够细（渐近区域），(3)解光滑。如果"2阶精确的"代码观测阶只有 0.8，基本肯定是边界条件用了1阶。这样 Richardson 外推就成了诊断工具——检验代码是否真的达到了宣称的精度。
🙋
默认值里 p=1.983 说是"2阶"。这是中点法则或 RK2 吗？
🎓
对，就是。这个例子用中点法则算 ∫₀¹ 1/(1+x²) dx = π/4 ≈ 0.7854。中点法则理论上2阶精确，观测结果 1.983 接近2。不过这个问题因为对称性产生了"超收敛"——中点法则实际会收敛到离真值略远的极限。这教会我们：观测阶只能告诉你收敛趋势和速度，u_ext 是否真的接近物理正确的答案还要另外验证。
🙋
在 CFD 里做"网格收敛"时，实际怎么用这个工具？
🎓
3步走：(1)制作3个网格——粗、中、细，节点数约成倍增加；(2)用同样物理和算法在每个网格上计算，记录代表物理量（阻力系数、应力峰值、努塞尔数等）为 u₃,u₂,u₁；(3)输入这个工具，得到 p 和 GCI。论文评审时，如果能证明"观测阶在理论阶的 ±0.5 以内，GCI 在允许范围"，就说明网格"已足够精细"。AIAA G-077 和 ASME V&V20 都把这个写成了标准。一个"只有漂亮图片"的单网格 CFD 报告在学术上是没有信用的。
🙋
渐近误差比 e₁/e₂ 应该接近"1/r^p"，偏离了怎么理解？
🎓
这是判断"是否进入渐近区域"的关键。如果 e₁/e₂ 与期望值（比如 r=2,p=2 时的 0.25）偏差大，说明3个点还不在渐近区域，p 的估计本身就不可靠。对策有两个：(a)取更细的网格，加第4个数据点做成成对比较；(b)检查解是否有不连续或强奇点。激波流、应力集中、重入角这些地方，解的光滑性会限制阶数，可能 p 永远卡在 1~1.5，这时候就"观测什么阶就报什么阶"，而不要硬推成理论阶。

常见问题

从3个步长 h₁<h₂<h₃（步长比 r=h₂/h₁=h₃/h₂）的数值解 u₁,u₂,u₃，用 p = ln(|u₂−u₃|/|u₁−u₂|) / ln(r) 计算。理论阶 p_theory（欧拉法=1、梯形法则或RK2=2、RK4=4）与观测值一致表示"代码和网格处于渐近区域"。观测阶与理论值偏差大，则怀疑：步长还是太粗、代码有 bug，或边界条件阶数较低。本工具即使步长比不等也按同一公式估算。
Richardson 外推值 u_ext = u₁ + (u₁−u₂)/(r^p − 1) 是消除渐近误差展开主项的推定值。3个解在"渐近区域"时，u_ext 的误差比原解的误差小 1~2 个数量级。判定标准是渐近误差比 e₁/e₂ 接近 1/r^p。例如 r=2, p=2 时期望 e₁/e₂ ≈ 0.25，若偏离较大则不应信任该阶数。Roache (1998) 的 GCI（网格收敛指数）用安全系数 Fs=1.25 定义误差置信带。
主要原因4个：(1)步长仍太粗，未入渐近区域——用更细步长重新计算。(2)边界条件/初始条件低阶实现（如1阶阶跃近似）——会支配整体离散误差，观测阶下降至0.5~1。(3)解在某处有奇点或不连续（激波、应力集中、凹角）——解的光滑性决定阶数，粗糙解达不到光滑解的理论阶。(4)迭代求解器不够收敛，迭代误差超过离散误差——需把残差下降至少1个数量级。
可以，正是此工具的主要用途。CFD 的 V&V 要求"用3网格做同样分析，报告观测阶、Richardson 外推值和 GCI"（AIAA G-077、ASME V&V20）。步骤：(1)粗、中、细3网格在同条件、同算法下计算；(2)代表物理量（阻力系数、峰值应力、Nu 数等）记为 u₃,u₂,u₁；(3)输入本工具确认 p、u_ext、GCI。若观测阶 p 与理论阶（FVM/FEM 通常2）的偏差在 ±0.5 内，判定为可信结果。

实际应用

CFD 的 V&V（验证与确认）：飞机阻力计算、涡轮机械效率、汽车外形空气动力等，商用 CFD（ANSYS Fluent、OpenFOAM、STAR-CCM+）的报告中观测收敛阶和 GCI 基本必须。AIAA G-077（航空航天）、ASME V&V20（流体热）、NAFEMS 指南都要求"至少3网格""确认观测阶""报告 GCI"。只有单网格"漂亮图"的报告在学术界会被否决。

有限元法的应力分析：应力集中部的峰值应力、J 积分等，网格细分后值波动明显。用粗、中、细3网格重新计算，用 Richardson 外推推定"真实应力"。重入角等奇点处，观测阶会跌到1以下，外推无法收敛，此时应"报告最细网格的值及其不确定性"。NAFEMS 基准问题多数包含这类收敛评估。

ODE 解法的精度评估：新时间积分方案（Runge-Kutta、Adams-Bashforth、IMEX 等）实现后，检验时间步 Δt 减半时误差是否按 2^p 下降。代码有 bug 时观测阶会是"奇怪的值"（本应3.2却是1.8），成为除错的强力线索。教科书习题和论文数值实验中频繁出现。

气象海洋模式：大气环流模式（GCM）和海洋模式（NEMO、MITgcm）也做空间/时间分辨率的收敛评估。但因强非线性和湍流，观测阶常低于理论值，进入渐近区域需庞大计算资源。实务上往往"收不了，但关键物理量不再变化"就停止，体现了 Richardson 外推的局限。

常见误区与注意事项

最常见的错误是"3点差值小就是收敛"的误解。比如 u₃=1.000, u₂=1.001, u₁=1.0011，u₃→u₂ 变 0.001，u₂→u₁ 仅 0.0001，计算出阶数 p=ln(10)/ln(2)≈3.3，然后欣喜"3阶精确！"，其实可能是"3点都粘在同一个错误的值（系统偏差）上了"。Richardson 外推是推定"收敛方向和速度"的工具，不保证"真值"。必须用另一方法（解析解、实验、其他求解器）交叉验证物理合理性。

其次，"不等步长仍用 p 公式"的错误。本工具的 p 公式默认 r=h₂/h₁=h₃/h₂ 相等。CFD 若用"自适应网格"或"混合细分"，步长比可能变 2/1.7，观测阶估计就会大幅偏差。实务中应用 Roache (1998) 的不等步长公式（迭代反解 p），或尽量做 r=2 的严格网格对。本工具只在 r≈const 时输出有意义的值。

第三，"阶数高就不用细网格"的误导。p=4（RK4 或高阶有限元）即便观测出来，现网格"误差够小"也是另一回事。误差形如 e≈C·h^p，系数 C 大的话即使 p 高，粗网格也会产生很大误差。GCI 公式中 E_h=(u₁−u₂)/(r^p−1) 是"当前实际的误差"，才是评判"够不够细"的标准，不能只看 p。"高阶手法在粗网格也高效"是进入渐近区域以后的事，初始时粗网格的低阶方法反而可能误差更小。

使用指南

  1. 输入粗、中、细3个步长 h₁、h₂、h₃（例：从 CFD 网格的单元数反推步长，推荐步长比 r=2）
  2. 输入各步长下的数值解 u₁、u₂、u₃（例：梁的最大挠度、压力降、热流通量）
  3. 模拟器自动计算观测收敛阶 p、Richardson 外推值 u_ext、截断误差 ε，判定数值方法的收敛性

具体计算示例

四点支撑梁的 FEM 分析，网格分割数 8→16→32 单元（步长比 r=2）时，最大挠度 u₁=10.24mm、u₂=10.12mm、u₃=10.08mm：p=(ln|u₁-u₂|/|u₂-u₃|)/ln(r)≈2.0，确认2阶精度。Richardson 外推值 u_ext≈10.07mm 为真值推定。RK2 法或 Crank-Nicolson 法可用同法验证阶数。

实务注意事项