Richardson 外推 | Codeboy's Blog

形式陈述

设 $A$ 是待求的精确量， $A(h)$ 是采用步长 $h>0$ 得到的近似。若其离散化误差在 $h\to0$ 时具有渐近展开

A(h)=A+ch^p+O(h^{p+q}), \qquad p,q>0,

其中 $c$ 与 $h$ 无关， $O(h^{p+q})$ 表示剩余项在足够小的 $h$ 下被某个常数乘 $h^{p+q}$ 控制，主误差随网格步长的幂次 $p$ 已由理论或独立验证确定。对步长 $h$ 与 $h/2$ ，组合

R(h)=\frac{2^pA(h/2)-A(h)}{2^p-1}

恰好消去 $ch^p$ ，从而得到

R(h)=A+O(h^{p+q}).

外推依赖的是相同极限、相同主误差系数和可靠的幂次 $p$ 。若两次计算改变了模型、边界处理、停止容差或随机样本，代数消项就没有共同误差展开可用。

当展开按 $p,p+q,p+2q,\ldots$ 递进，取几何步长 $h_j=h_0/2^j$ ，可构造外推表

R_{j,0}=A(h_j),

R_{j,k} =R_{j,k-1} +\frac{R_{j,k-1}-R_{j-1,k-1}} {2^{p+(k-1)q}-1}, \qquad 1\le k\le j.

每升一列消去一个已知幂次。生成深度 $m$ 的完整表需要 $m+1$ 次基础近似；表内组合为 $O(m^2)$ 次标量运算，若逐行计算并用一维数组覆盖上一行，则存储可降为 $O(m)$ 。总成本通常由最细步长的 $A(h_m)$ 决定，必须与它的网格或函数求值成本一起报告。

在单主项模型下，细网格近似自身的误差可估为

|A(h/2)-A| \approx \frac{|A(h/2)-A(h)|}{2^p-1}.

这是渐近误差估计，不是无条件严格上界。实际停止可要求最新对角线变化小于绝对—相对混合容差，同时检查相邻差值比接近预期的 $2^p$ ；还应设置最大层数，并在差值增大、符号无规律翻转或步长已到浮点平台时退出失败。

直觉

外推把误差当作可以辨认的信号。若粗、细两个结果都带有同一种 $h^p$ 偏差，只是细网格把它缩小了 $2^p$ 倍，就能选择一组权重让这份偏差相消，留下下一阶项。它不是从两个相近数字中凭空制造精度，而是兑现一条已经知道的误差展开。

进入渐近区之前，多个误差项可能同量级；进入舍入平台之后，差值又主要反映数值噪声。只有两者之间的区间，按理论比例缩放的主误差项才清晰可辨。

例子与边界

对足够光滑的函数，中心差分 $D(h)=[f(x+h)-f(x-h)]/(2h)$ 满足

D(h)=f'(x)+c_2h^2+c_4h^4+O(h^6).

取 $p=q=2$ ，第一次外推得到

R(h)=\frac{4D(h/2)-D(h)}{3} =f'(x)+O(h^4).

可以完整算一个多项式例子。取 $f(t)=t+t^3+t^5$ ，在 $x=0$ 处真导数为 $1$ ，而 $D(h)=1+h^2+h^4$ 。因此 $D(1/2)=21/16$ ， $D(1/4)=273/256$ ，外推得到 $R(1/2)=63/64$ 。细步长原误差为 $17/256$ ，外推误差变成 $-1/64$ ；代数上 $R(h)=1-h^4/4$ ，所以继续把 $h$ 减半会让绝对误差缩小 $16$ 倍，而非原方法的渐近 $4$ 倍。

这里的四阶不是由两次结果“看起来接近”得出，而是由中心差分只含偶次主误差的展开保证。若边界 stencil 混入 $O(h^3)$ 项，仍套用偶次外推表便会给出错误阶数。

若误差形如 $h^p\log h$ ，用常数比例 $2^p$ 不能完整消去主项；若误差随 $1/h$ 振荡，相邻网格甚至可能偶然相等。函数不够光滑、网格跨过间断或迭代求解误差大于离散误差时，也不会出现所需的稳定比值。

外推不能自动消除基础计算中的噪声。若 $A(h)$ 与 $A(h/2)$ 各有绝对误差至多 $\varepsilon$ ，一次组合带来的噪声至多为 $(2^p+1)\varepsilon/(2^p-1)$ ； $p=2$ 时是 $5\varepsilon/3$ 。系数并非必然随层数无限增大，真正需要检查的是主截断误差已被消去后，舍入或求解噪声是否开始占主导。继续填表到对角线“数字不再变化”可能只是舍入停滞；当误差估计停止按预期阶数下降时，应终止而不是无限增加层数。

推论与应用

这里的 $p$ 是 $h\to0$ 的网格误差幂次，不能当作迭代 Q-阶。例如 $A(h)=A+h^2$ ，取 $h_j=2^{-j}$ ，误差 $e_j=4^{-j}$ 满足 $e_{j+1}/e_j=1/4$ ，所以这列近似是 Q-线性，虽然网格误差为二阶。离散化误差说明 $A(h)$ 的误差来源。Romberg 积分把复合梯形规则的偶次幂展开代入本页递推，ODE 步长加倍也使用相同消项结构；每个应用仍需单独证明自己的展开和稳定区间。

可审计的外推结果应保存原始 $A(h_j)$ 、假定的 $p,q$ 、外推表、差值比和退出原因。只报告最后一格数字，会掩盖方法是否真正处在渐近区。

参考资料

Robert Israel, “Richardson Extrapolation”, University of British Columbia course notes, 2002.
Lewis Fry Richardson, “The Approximate Arithmetical Solution by Finite Differences of Physical Problems,” Philosophical Transactions of the Royal Society A 210, 1911.
Josef Stoer and Roland Bulirsch, Introduction to Numerical Analysis, 3rd ed., Springer, 2002, extrapolation methods.

形式陈述 ​

反向链接

形式陈述