理查森外推

理查森外推（Richardson extrapolation）组合不同步长得到的近似值，消去一个已知阶数的主导误差项。它依赖误差展开中稳定的步长幂次，适用于数值微分、积分、常微分方程及其他网格近似；仅知道“步长较小通常更准”还不足以决定组合系数。

例如，用中心差分计算 $f(x)=x+x^3$ 在零处的导数。真值为一，而步长 $h$ 的差商为 $$D(h)=\frac{f(h)-f(-h)}{2h}=1+h^2.$$ 取 $h=1$ 得二，$h=1/2$ 得 $5/4$。二次误差在步长减半后变为四分之一，因此组合 $(4D(1/2)-D(1))/3$ 恰好消去它，得到一。

从误差模型推导系数

设某个近似满足 $$A(h)=A+c{h}^p+O(h^{p+q}),p>0,q>0,$$ 其中 $A$ 为目标值，$c$ 不随 $h$ 变化。固定缩放比 $r>1$，则 $$A(h/r)=A+c{r}^{-p}{h}^p+O(h^{p+q}).$$ 将细步长项乘 $r^p$ 再减去粗步长项，$c{h}^p$ 正好抵消；常数项剩 $(r^p-1)A$。于是 $$R(h)=\frac{r^{p}A(h/r)-A(h)}{r^p-1}=A+O(h^{p+q}).$$ 这就是外推公式的完整推导。它没有要求知道 $c$ 的值，却要求阶数 $p$ 正确，且两次近似处在同一误差展开有效的范围内。

同一模型也说明细网格误差近似为 $$A(h/r)-A\approx \frac{A(h)-A(h/r)}{r^p-1}.$$ 外推值相对细网格值的改正量恰为这个估计的相反数。如果主误差系数恰好为零，或误差被舍入主导，这个近似关系不能提供预期的可靠估计。

中心差分的二次主误差

若 $f$ 在查询点附近足够光滑，泰勒展开为 $$\begin{array}{cc}f(x+h)& =f(x)+h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{″}(x)+\frac{h^3}{6}{f}^{‴}(x)+O(h^4),\\ f(x-h)& =f(x)-h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{″}(x)-\frac{h^3}{6}{f}^{‴}(x)+O(h^4).\end{array}$$ 两式相减再除 $2h$，首先得到二次主误差；若继续展开到五阶并控制余项，还可得到 $D(h)=f^{\prime }(x)+h^2{f}^{‴}(x)/6+O(h^4)$。所以 $p=2,r=2$ 时系数为四，外推得到四阶近似，而不是套用前向差分常用的 $p=1$。

主例误差恰与h平方成正比，外推相当于取直线在零步长处的截距。

图中横坐标为 $h^2$，并非 $h$。两次取值相当于点 $(1,2)$ 与 $(1/4,5/4)$；直线截距是一，和代数消元一致。

阶数错误导致变差的变式

变式。 将函数换成 $f(x)=x+x^5$，仍在零处取中心差分。这次 $D(h)=1+h^4$，二次系数恰为零；两个步长结果为 $D(1)=2,D(1/2)=17/16$。若仍用 $p=2$，外推值为 $$\frac{4(17/16)-2}{3}=\frac{3}{4},$$ 误差四分之一，反比细差分的误差十六分之一更大。若依据这里已知的精确误差使用 $p=4$，则 $(16(17/16)-2)/15=1$。在一般数值问题中，不能因两次结果接近，就断定选用的阶数一定正确。

需要估计观测阶时，至少可比较三个连续步长：若主项确实非零且占主导，$|A(h)-A(h/r)|/|A(h/r)-A(h/r^2)|$ 应接近 $r^p$。主例差值为 $3/4$ 与 $3/16$，比值四，支持二阶模型；差值为零时这一检验没有定义，不能强行取对数。

可运行的组合与中心差分

组合函数返回外推值，以及按主误差模型估计的细步长误差绝对值。它不证明误差模型成立。所有输入必须有限，步长比大于一、阶数为正；输入不被改写。

Python 3

from math import isfinite

def richardson(coarse, fine, ratio, order):
    coarse, fine, ratio, order = map(float, (coarse, fine, ratio, order))
    if not all(isfinite(z) for z in (coarse, fine, ratio, order)):
        raise ValueError("finite input required")
    if ratio <= 1 or order <= 0:
        raise ValueError("ratio must exceed one and order must be positive")
    try:
        factor = ratio**order
    except OverflowError as error:
        raise ArithmeticError("extrapolation factor overflow") from error
    if not isfinite(factor) or factor == 1:
        raise ArithmeticError("unrepresentable extrapolation factor")
    result = (factor*fine-coarse)/(factor-1)
    estimate = abs(fine-coarse)/(factor-1)
    if not isfinite(result) or not isfinite(estimate):
        raise ArithmeticError("nonfinite extrapolated result")
    return result, estimate

def central_difference(function, x, h):
    x, h = float(x), float(h)
    if not isfinite(x) or not isfinite(h) or h <= 0:
        raise ValueError("finite point and positive step required")
    if x+h == x or x-h == x:
        raise ArithmeticError("step lost to floating-point rounding")
    upper, lower = float(function(x+h)), float(function(x-h))
    if not isfinite(upper) or not isfinite(lower):
        raise ValueError("nonfinite function value")
    return (upper-lower)/(2*h)

f = lambda x: x+x**3
coarse, fine = central_difference(f, 0, 1), central_difference(f, 0, 0.5)
assert (coarse, fine) == (2, 5/4)
assert richardson(coarse, fine, 2, 2) == (1, 1/4)
g = lambda x: x+x**5
coarse, fine = central_difference(g, 0, 1), central_difference(g, 0, 0.5)
assert richardson(coarse, fine, 2, 2)[0] == 3/4
assert richardson(coarse, fine, 2, 4)[0] == 1
print(richardson(2, 5/4, 2, 2))

已有两个近似时，外推组合需 $O(1)$ 时间和空间；中心差分每个步长需要两次函数评价。若近似来自大型网格计算，主要成本仍在生成两次近似，不能因为组合只需几次算术就称整个计算为常数成本。

舍入与光滑性边界

若两个输入近似各含绝对误差至多 $\epsilon$，外推中的线性组合可能把这部分误差放大至 $(r^p+1)\epsilon /(r^p-1)$。例如二阶减半外推的系数绝对值之和为 $5/3$。消去截断误差并不同时消去舍入噪声，尤其步长过小时差分中的相消会更严重。

非光滑函数也可能给出误导性的稳定数字。对 $|x|$ 在零处做中心差分，任意步长都得到零，外推也得到零，但该点导数本来就不存在。外推不能补回近似问题中缺失的数学条件。