Richardson 外推法

概述

Richardson 外推法（Richardson Extrapolation）由刘易斯·弗赖·理查森于 1911 年发表（1927 年系统阐述为"延迟趋近极限法"）。其核心思想是：数值方法的误差通常具有渐近展开式 $A(h) = A^ + a_1 h^p + O(h^{p+1})$，通过组合两个不同步长下的计算结果，可以代数消去主导误差项 $a_1 h^p$*，将精度从 $O(h^p)$ 提升到 $O(h^{p+1})$，同时无偿获得误差估计。外推法是数值分析中最通用的"元方法"——可叠加在几乎任何具有渐近误差展开的数值方法之上。

关键内容

核心公式

假设：误差渐近展开 $A(h) = A^* + a_1 h^p + O(h^{p+1})$（$p$ 是方法阶数，$a_1$ 与 $h$ 无关）

步长加倍外推（$r=2$，最常用）：

$$A_{\text{improved}} = \frac{2^p \cdot A(h/2) - A(h)}{2^p - 1}$$

误差从 $O(h^p)$ 提升到 $O(h^{p+1})$——精度提高一阶。

一般步长比 $r$：

$$A_{\text{improved}} = \frac{r^p \cdot A(h/r) - A(h)}{r^p - 1}$$

代数本质：类比方程组消元法——$A(h)$ 和 $A(h/2)$ 是两个含 $A^$ 和 $a_1$ 的"方程"，通过适当线性组合消去 $a_1$（无需知道其数值），直接得到对 $A^$ 的更优估计。

递归外推表

外推可以递归应用：先消去 $O(h^p)$ 项，再消去 $O(h^{p+1})$ 项……每层提高一阶精度。

设 $T_0^{(k)} = A(h/2^k)$，递推：

$$T_j^{(k)} = T_{j-1}^{(k+1)} + \frac{T_{j-1}^{(k+1)} - T_{j-1}^{(k)}}{2^{jp} - 1}$$

$T_0^{(0)}$
$T_0^{(1)}$	$T_1^{(0)}$
$T_0^{(2)}$	$T_1^{(1)}$	$T_2^{(0)}$

对角线方向精度依次提升一阶。这就是 Romberg积分的核心结构。

误差估计：免费副产品

$$A(h) - A^* \approx \frac{A(h/2) - A(h)}{2^p - 1} \cdot 2^p$$

这给出了对主导误差的量化估计，是现代自适应数值方法中最常用的误差指示器之一。

数字示例（梯形法则积分 $\int_0^1 e^x dx = e-1$，$p=2$）

步长 $h$	梯形法则误差	外推后误差
1.0 → 0.5	$0.141$ → $0.036$	$0.00059$（降低约 60 倍）
0.5 → 0.25	$0.036$ → $0.0089$	$0.000034$

一次外推将误差阶从 $O(h^2)$ 提升至 $O(h^4)$。

适用条件与局限性

需要误差渐近展开成立：当 $h$ 太大（高阶项不可忽略）或解有奇异性时，标准展开可能不成立。

舍入误差限制：两结果相减可能发生数值抵消（cancellation），实际精度通常不超过机器精度的平方根。

计算成本：需要多个步长下的计算；高维 PDE 问题步长减半则计算量 $\sim 2^d$ 增长，需权衡收益。

不稳定格式无效：无论多好的误差分析都无法挽救数值不稳定的格式（理查森自己的天气预报失败正是这一教训）。

主要后续发展

方法	年份	关系
Romberg积分	1955	递归外推应用于梯形法则，利用欧拉-麦克劳林展开只有偶数次幂的特殊性
Bulirsch-Stoer 方法	1966	外推法用于 ODE 初值问题，天体力学精度首选
嵌入式 Runge-Kutta	1960s+	步长控制基于本质上等同于 Richardson 误差估计的嵌入式对
CFD 网格收敛指数（GCI）	1994	外推法的标准化 CFD 验证工具（Roache 1994/1998）
自适应网格加密（AMR）	1980s+	基于 Richardson 误差估计驱动的自适应网格精化

来源

• raw/books/数值分析/10_richardson_extrapolation.md

相关

• 刘易斯·弗赖·理查森
• Romberg积分
• 高斯求积公式
• 欧拉方法
• Jacobi迭代法