龙贝格积分的外推加速原理及Richardson外推法的一般框架

计算积分时，我们常用梯形法则。它的想法很简单：把一条曲线下的面积，近似成一连串小梯形面积的总和。分的梯形越多，每个梯形越窄，结果就越接近真实值。这个过程可以不断二分区间，得到一个梯形序列 $T_1, T_2, T_3, ...$，每个都比前一个更精确。

核心问题：我们能否利用这个已有的序列，聪明地组合它们，从而直接得到一个比其中任何一项都更精确的近似值？

答案是肯定的。这就是龙贝格积分（Romberg Integration）和Richardson外推法的精髓。

第一部分：理解梯形法则的误差

要“聪明地组合”，首先要明白我们组合的对象存在什么“误差”。

梯形法则的误差，主要来自用直的梯形边去代替弯曲的函数曲线。当二分区间的次数为 $k$（即区间宽度为 $h = \frac{b-a}{2^k}$）时，梯形法则结果 $T_k$ 的误差可以展开为关于步长 $h$ 的幂级数：

$$T_k = I + a_1 h^2 + a_2 h^4 + a_3 h^6 + ...$$

其中 $I$ 是我们想要求的精确积分值。$a_1, a_2, ...$ 是与被积函数有关的常数，但与 $k$ 无关。

关键发现：误差的主要项是 $a_1 h^2$。这意味着，当我们将步长 $h$ 减半（从 $h$ 变为 $h/2$，即从 $T_k$ 变为 $T_{k+1}$），新的误差项会变为 $a_1 (h/2)^2 = \frac{1}{4} a_1 h^2$。

这表明，简单的梯形序列虽然收敛，但收敛速度是“2阶”的（误差随 $h^2$ 减小）。我们能否消除 $a_1 h^2$ 这个主要误差项，从而获得一个误差为 $O(h^4)$ 甚至更高阶的新序列？这便是外推加速的目标。

第二部分：Richardson外推法的一般框架

Richardson外推法是一个通用框架，适用于任何误差具有已知幂级数形式的计算方法。它通过组合不同精度（通常对应不同步长）的计算结果，系统性地消除误差展开式中的低阶项。

通用公式推导：

假设一个计算方法 $M(h)$ 在步长为 $h$ 时对真值 $M^*$ 的近似满足：
$$M(h) = M^* + c_1 h^{p_1} + c_2 h^{p_2} + c_3 h^{p_3} + ...$$
其中 $0 < p_1 < p_2 < ...$。

目标：找到一个组合 $M(h)$ 和 $M(h/t)$ 的新方法，消除主要误差项 $c_1 h^{p_1}$。这里 $t$ 是一个大于1的常数（常取 $t=2$，即步长减半）。

执行推导：

1. 写出两个误差表达式：
   $$M(h) = M^* + c_1 h^{p_1} + c_2 h^{p_2} + ...$$
   $$M(h/t) = M^* + c_1 (h/t)^{p_1} + c_2 (h/t)^{p_2} + ...$$

2. 消除 $c_1 h^{p_1}$ 项。将第一个式子乘以 $t^{p_1}$，与第二个式子相减：
   $$t^{p_1} M(h) - M(h/t) = (t^{p_1}-1)M^* + [t^{p_1}c_2 h^{p_2} - c_2 (h/t)^{p_2}] + ...$$

3. 解出 $M^*$ 的新近似 $M_{new}$：
   $$M_{new} = \frac{t^{p_1} M(h) - M(h/t)}{t^{p_1} - 1}$$

结论：这个新的 $M_{new}$ 的误差主项已经变为 $O(h^{p_2})$，比原来的 $O(h^{p_1})$ 更高阶。这意味着我们一次外推，就将收敛速度从 $p_1$ 阶提升到了 $p_2$ 阶。

第三部分：将Richardson外推应用于梯形法则——龙贝格积分

将Richardson外推的通用框架应用到梯形法则上，就得到了龙贝格积分。

具体应用：

1. 确认误差形式：对于梯形法则 $T(h)$，误差主项是 $h^2$ 项。即 $p_1 = 2$， $t=2$。

2. 套用Richardson公式：用 $T(h)$ 和 $T(h/2)$ 组合，消除 $h^2$ 项。
   $$R_{new} = \frac{2^{2} \cdot T(h) - T(h/2)}{2^{2} - 1} = \frac{4 T(h) - T(h/2)}{3}$$

3. 解释新结果：这个 $R_{new}$ 被称为辛普森法则的近似值 $S(h)$。它的误差主项已经是 $O(h^4)$ 了。我们从梯形序列 $T_k$ 出发，通过一次简单的加权平均，直接得到了一个精度高得多的辛普森序列 $S_k$。

4. 继续外推：既然 $S(h)$ 的误差主项是 $h^4$ 项（即 $p_1=4$），我们能否对辛普森序列再次应用Richardson外推，消除 $h^4$ 项？完全可以。
   用 $S(h)$ 和 $S(h/2)$ 组合（此时 $p_1=4$）：
   $$B_{new} = \frac{2^{4} \cdot S(h) - S(h/2)}{2^{4} - 1} = \frac{16 S(h) - S(h/2)}{15}$$
   这个 $B_{new}$ 就是布尔法则的近似值 $B(h)$，误差主项是 $O(h^6)$。

龙贝格积分表格的构建：

上述过程可以系统性地进行，形成一个表格（称为龙贝格表格）。表格的列代表外推的阶数，行代表二分区间的次数。

构建步骤：

1. 定义梯形序列。

   • 第一列 ($R_{i,1}$) 存储梯形法则的结果 $T_k$。例如，$R_{1,1} = T_1$ (用1个梯形), $R_{2,1} = T_2$ (用2个梯形), ...。

2. 应用Richardson外推递推公式。

   • 后续列的每一项 $R_{i,j}$ ($j > 1$) 都由前一列的两项 $R_{i,j-1}$ 和 $R_{i-1,j-1}$ 计算得出。
   • 通用递推公式：
     $$R_{i,j} = R_{i,j-1} + \frac{R_{i,j-1} - R_{i-1,j-1}}{4^{j-1} - 1}$$
     这个公式的本质就是Richardson外推。分母 $4^{j-1}$ 来源于 $t^{p_1}$，其中 $t=2$，而当前要消除的误差项阶数为 $p_1 = 2 \times (j-1)$。

3. 观察表格对角线。

   • 表格对角线上的元素 $R_{i,i}$ 代表不断外推后得到的最佳估计。$R_{1,1}$ 是梯形法则，$R_{2,2}$ 是辛普森法则，$R_{3,3}$ 是布尔法则，以此类推。沿着对角线向下移动，近似值的精度会阶跃式提高。

第四部分：操作指南与示例

假设我们要计算积分 $I = \int_{0}^{1} e^x dx$ 的近似值（真值为 $e - 1 \approx 1.718281828$）。

第一步：生成梯形序列 ($R_{i,1}$)。

1. 计算 $T_1$ ($R_{1,1}$)：区间 [0,1]， 1个梯形。
   $$T_1 = \frac{1-0}{2} [f(0) + f(1)] = 0.5 \times (e^0 + e^1) \approx 0.5 \times (1 + 2.71828) = 1.85914$$

2. 计算 $T_2$ ($R_{2,1}$)：将区间二分，2个梯形。
   中点 $x=0.5$， $f(0.5)=e^{0.5}\approx 1.64872$。
   $$T_2 = \frac{0.5}{2}[f(0) + 2f(0.5) + f(1)] \approx 0.25 \times (1 + 3.29744 + 2.71828) \approx 1.75393$$

3. 计算 $T_3$ ($R_{3,1}$)：再二分，4个梯形。
   新增点 $x=0.25$ ($f \approx 1.28403$), $x=0.75$ ($f \approx 2.11700$)。
   $$T_3 \approx 1.71886$$

第二步：应用Richardson外推，填充表格的下一列 ($R_{i,2}$)。

1. 计算 $R_{2,2}$ (辛普森法则 $S_1$)：
   $$R_{2,2} = R_{2,1} + \frac{R_{2,1} - R_{1,1}}{4^{1} - 1} = 1.75393 + \frac{1.75393 - 1.85914}{3} \approx 1.75393 - 0.03507 = 1.71886$$

2. 计算 $R_{3,2}$ (辛普森法则 $S_2$)：
   $$R_{3,2} = R_{3,1} + \frac{R_{3,1} - R_{2,1}}{3} \approx 1.71886 + \frac{1.71886 - 1.75393}{3} \approx 1.71886 - 0.01169 = 1.71717$$

第三步：继续外推，填充更高级列 ($R_{i,3}$ 等)。

1. 计算 $R_{3,3}$ (布尔法则 $B_1$)：
   $$R_{3,3} = R_{3,2} + \frac{R_{3,2} - R_{2,2}}{4^{2} - 1} = 1.71717 + \frac{1.71717 - 1.71886}{15} \approx 1.71717 - 0.00011 = 1.71706$$

第四步：检查收敛。

• 真值： $1.718281828$
• $R_{1,1}$ (梯形)： $1.85914$ (误差约 8%)
• $R_{2,2}$ (辛普森)： $1.71886$ (误差约 0.03%)
• $R_{3,3}$ (布尔)： $1.71706$ (误差约 0.07%)

在这个小例子中，$R_{3,3}$ 略逊于 $R_{2,2}$，是因为我们使用的数据点太少（只有4个），高阶方法的优势无法充分发挥。在实际应用中，随着数据点的增加，龙贝格对角线上的高阶项会迅速收敛到远高于低阶项的精度。

总结龙贝格积分的执行流程：

1. 选择被积函数和积分区间。
2. 计算初始的梯形序列 $T_1, T_2, T_4, ...$（即 $R_{i,1}$）。通常通过不断二分区间并利用已有计算结果来高效完成。
3. 设置一个空的二维表格，第一列为 $R_{i,1}$。
4. 循环：对于每一列 $j$ (从2开始)，使用上一列的数据 $R_{i,j-1}$ 和 $R_{i-1,j-1}$，通过递推公式 $R_{i,j} = R_{i,j-1} + (R_{i,j-1} - R_{i-1,j-1}) / (4^{j-1} - 1)$ 计算当前列的所有元素 $R_{i,j}$。
5. 重复步骤4，直到计算到所需的列数或直到表格对角线上的相邻值之差小于预设的容差 $\epsilon$。
6. 输出表格对角线或最后一列中满足精度要求的值作为最终积分近似。