Richardson外推

Richardson外推（Richardson Extrapolation）是一种威力强大的数值技巧，它能把低精度结果”加工”成高精度结果。核心思想是利用截断误差的渐近展开，通过组合不同步长的计算结果来消去误差主导项。

---

一、最直观的例子：计算导数

问题：计算 $f^{\prime }(x)$，已知 $f(x)$

方法1：向前差商（一阶）

$$D_{h}f(x)=\frac{f(x+h)-f(x)}{h}$$

误差分析（泰勒展开）：

$$f(x+h)=f(x)+h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{\prime \prime }(x)+\frac{h^3}{6}{f}^{\prime \prime \prime }(x)+\cdots$$

代入：

$$D_{h}f(x)=f^{\prime }(x)+\frac{h}{2}{f}^{\prime \prime }(x)+\frac{h^2}{6}{f}^{\prime \prime \prime }(x)+\cdots$$

误差 = $O(h)$

方法2：中心差商（二阶）

$$C_{h}f(x)=\frac{f(x+h)-f(x-h)}{2h}$$

泰勒展开：

$$f(x+h)=f(x)+h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{\prime \prime }(x)+\frac{h^3}{6}{f}^{\prime \prime \prime }(x)+\frac{h^4}{24}{f}^{(4)}(x)+\cdots$$ $$f(x-h)=f(x)-h{f}^{\prime }(x)+\frac{h^2}{2}{f}^{\prime \prime }(x)-\frac{h^3}{6}{f}^{\prime \prime \prime }(x)+\frac{h^4}{24}{f}^{(4)}(x)-\cdots$$

相减：

$$C_{h}f(x)=f^{\prime }(x)+\frac{h^2}{6}{f}^{\prime \prime \prime }(x)+O(h^4)$$

误差 = $O(h^2)$（更好！）

但是：中心差商需要 $f(x-h)$，有时不方便。

---

二、Richardson的洞见：用两个低阶结果构造高阶结果

假设我们有一个近似 $A(h)$，它逼近真值 $A^{\ast }$：

$$A(h)=A^{\ast }+C_{1}h+C_2{h}^2+C_3{h}^3+\cdots$$

其中 $C_k$ 是常数（与 $h$ 无关）。

关键步骤：

1. 用步长 $h$ 计算：$A(h)=A^{\ast }+C_{1}h+C_2{h}^2+\cdots$
2. 用步长 $h/2$ 计算：$A(h/2)=A^{\ast }+C_1(h/2)+C_2(h/2{)}^2+\cdots$

消去 $h$ 项：

观察：

$$2A(h/2)=2A^{\ast }+C_{1}h+C_2\frac{h^2}{2}+\cdots$$ $$A(h)=A^{\ast }+C_{1}h+C_2{h}^2+\cdots$$

相减：

$$2A(h/2)-A(h)=(2A^{\ast }-A^{\ast })+(C_{1}h-C_{1}h)+\left(\frac{C_2{h}^2}{2}-C_2{h}^2\right)+\cdots$$ $$=A^{\ast }-\frac{C_2{h}^2}{2}+\cdots$$

神奇的事情发生了：

• 新的组合 $2A(h/2)-A(h)$ 的误差是 $O(h^2)$
• 而原来的 $A(h)$ 误差是 $O(h)$
• 我们只用加法和乘法，就把精度提高了一阶！

定义外推结果：

$$A_1(h)=2A(h/2)-A(h)$$

现在 $A_1(h)=A^{\ast }+O(h^2)$

---

三、一般公式：Richardson外推表

假设近似满足更一般的渐近展开：

$$A(h)=A^{\ast }+C_1{h}^p+C_2{h}^{p+1}+C_3{h}^{p+2}+\cdots$$

其中 $p$ 是方法的阶。

外推公式：

$$A_1(h)=\frac{2^{p}A(h/2)-A(h)}{2^p-1}$$

验证：

$$A(h)=A^{\ast }+C_1{h}^p+O(h^{p+1})$$ $$A(h/2)=A^{\ast }+C_1(h/2{)}^p+O(h^{p+1})$$ $$2^{p}A(h/2)=2^p{A}^{\ast }+C_1{h}^p+O(h^{p+1})$$ $$2^{p}A(h/2)-A(h)=(2^p-1)A^{\ast }+O(h^{p+1})$$ $$A_1(h)=A^{\ast }+O(h^{p+1})$$

---

四、具体数值例子

例子1：计算 $\pi$ 的经典方法（圆周率）

用圆面积公式：$\pi =4\times$（单位圆第一象限面积）

数值积分近似：把 $[0,1]$ 分成 $n$ 份，用矩形法：

方法：$A(n)=\frac{4}{n}{\sum }_{i=1}^n\sqrt{1-{\left(\frac{i-0.5}{n}\right)}^2}$

• $n=4$（大步长 $h=1/4$）：

  $$A(4)=1\times [\sqrt{1-(0.125{)}^2}+\sqrt{1-(0.375{)}^2}+\sqrt{1-(0.625{)}^2}+\sqrt{1-(0.875{)}^2}]$$ $$\approx 1\times [0.9922+0.9270+0.7806+0.4841]=3.1839$$

• $n=8$（步长减半 $h=1/8$）：

  $$A(8)\approx 3.1395$$

Richardson外推（矩形法阶数 $p=1$？等等…）

实际上矩形法对光滑被积函数是二阶的（中点矩形公式），所以 $p=2$：

$$\pi \approx \frac{4\times 3.1395-3.1839}{4-1}=\frac{12.5580-3.1839}{3}=\frac{9.3741}{3}=3.1247$$

还不够好，但比原来好。

其实：中点矩形公式误差展开：

$$A(h)=A^{\ast }+C_1{h}^2+C_2{h}^4+\cdots$$

所以 $p=2$，外推后误差 $O(h^4)$。

---

例子2：欧拉法解 ODE 的外推

解 $y^{\prime }=-y,y(0)=1$，求 $y(1)\approx e^{-1}=0.36787944$

显式欧拉法（$p=1$）：

步长 $h=1$（一步到位）：

$$y_1=1+1\times (-1)=0$$

步长 $h=0.5$（两步）：

• 第一步：$y_{0.5}=1+0.5\times (-1)=0.5$
• 第二步：$y_1=0.5+0.5\times (-0.5)=0.25$

Richardson外推（$p=1$）：

$$y_{\text{extrap}}=\frac{2^1\times 0.25-0}{2^1-1}=0.5-0=0.5$$

还是不准，但比 0 和 0.25 都好。

再细化：
步长 $h=0.25$（四步）：

• 计算得 $y_1\approx 0.31640625$

用 $h=0.5$ 和 $h=0.25$ 外推：

$$y_{\text{extrap}}=\frac{2\times 0.31640625-0.25}{1}=0.6328125-0.25=0.3828125$$

误差：$0.3828125-0.367879=0.014933$

再与更小步长组合可以继续改进。

---

五、递推外推：Romberg积分法

最著名的应用是Romberg积分法——梯形法的Richardson外推。

梯形公式误差展开：

$$T(h)=I^{\ast }+C_1{h}^2+C_2{h}^4+C_3{h}^6+\cdots$$

其中 $I^{\ast }$ 是精确积分值。

Romberg表构造：

设 $T_{k,0}$ 是将区间 $2^k$ 等分的梯形公式结果。

外推表：

$$\begin{array}{ccccc}h& T_{0,0}& & & \\ h/2& T_{1,0}& T_{1,1}& & \\ h/4& T_{2,0}& T_{2,1}& T_{2,2}& \\ h/8& T_{3,0}& T_{3,1}& T_{3,2}& T_{3,3}\end{array}$$

其中：

$$T_{k,m}=\frac{4^m{T}_{k,m-1}-T_{k-1,m-1}}{4^m-1},m=1,2,\dots ,k$$

解释：

• $T_{k,0}$：梯形公式，误差 $O(h^2)$
• $T_{k,1}$：一次外推，误差 $O(h^4)$
• $T_{k,2}$：二次外推，误差 $O(h^6)$
• 以此类推

例子：计算 ${\int }_0^1\frac{4}{1+x^2}dx=\pi$

k	$T_{k,0}$（梯形）	$T_{k,1}$（Simpson）	$T_{k,2}$	$T_{k,3}$
0	3.00000000
1	3.10000000	3.13333333
2	3.13117647	3.14156863	3.14211765
3	3.13898849	3.14159250	3.14159409	3.14158578
4	3.14094161	3.14159265	3.14159266	3.14159264

注意：$T_{4,3}$ 已经达到 $3.14159264$，与 $\pi$ 相差 $10^{-8}$ 量级！

---

六、Richardson外推的数学原理

渐近展开假设

Richardson外推成立的关键是假设近似值可以展开为：

$$A(h)=A^{\ast }+{\alpha }_1{h}^{p_1}+{\alpha }_2{h}^{p_2}+{\alpha }_3{h}^{p_3}+\cdots$$

其中 $0<p_1<p_2<p_3<\cdots$，且 ${\alpha }_i$ 与 $h$ 无关。

一般外推公式

假设我们有两个步长 $h$ 和 $h/t$（通常 $t=2$）：

$$\{\begin{array}{l}A(h)=A^{\ast }+{\alpha }_1{h}^{p_1}+{\alpha }_2{h}^{p_2}+\cdots \\ A(h/t)=A^{\ast }+{\alpha }_1(h/t{)}^{p_1}+{\alpha }_2(h/t{)}^{p_2}+\cdots \end{array}$$

消去 ${\alpha }_1{h}^{p_1}$ 项：

$$A^{\ast }(h)=\frac{t^{p_1}A(h/t)-A(h)}{t^{p_1}-1}=A^{\ast }+{\beta }_2{h}^{p_2}+\cdots$$

递推外推（Richardson表）

我们可以重复这一过程，消除越来越多的误差项。

设 $A_{i,0}=A(h/t^i)$，则：

$$A_{i,j}=\frac{t^{p_j}{A}_{i,j-1}-A_{i-1,j-1}}{t^{p_j}-1},j=1,2,\dots ,i$$

其中 $A_{i,j}$ 的误差阶为 $h^{p_{j+1}}$。

---

七、实际应用的注意事项

1. 渐近展开必须存在

不是所有数值近似都有这样的展开。需要：

• 解足够光滑
• 离散格式规则

2. 舍入误差限制

外推不能无限进行，因为：

• 步长太小时，舍入误差主导
• 外推会放大舍入误差

经验法则：外推通常进行 3-5 次效果最好。

3. 计算成本 vs 精度增益

外推需要多次计算（不同步长），但通常比单纯减小步长更高效。

4. 自动检测展开阶数

有时阶数 $p$ 未知，可以通过三个不同步长的结果估计。

八、在ODE求解中的应用：Gragg-Bulirsch-Stoer方法

这是Richardson外推在ODE中的高级应用。

基本思想：

1. 用中点法（leapfrog）等偶数阶方法积分固定步数
2. 用不同步长（不同外推点）计算
3. 对结果进行多项式或有理函数外推
4. 得到高精度结果

优点：可以达到非常高的精度阶（如16阶以上）。

---

九、总结：Richardson外推的威力

核心公式（记住这个！）：

$$\text{高阶近似}=\frac{2^p\times (\text{步长减半结果})-(\text{原结果})}{2^p-1}$$

三个关键点：

1. 利用误差结构：假设误差有系统性的 $h^p$ 展开
2. 消去主导误差：通过线性组合消去最低阶误差项
3. 精度免费提升：只做算术运算，不增加函数求值次数（相对于已经算好的不同步长结果）

应用领域：

1. 数值积分：Romberg积分法
2. 数值微分：提高差分公式精度
3. ODE/PDE求解：步长减半误差估计和外推
4. 级数加速：收敛加速技术
5. 任何有渐近误差展开的数值计算

哲学意义：Richardson外推体现了数值分析的精髓——理解误差的结构，然后用智慧（而非蛮力）消除它。就像用两张模糊照片合成一张清晰照片，利用了模糊的规律性。