复合梯形公式的 Richardson 外推加速收敛技巧
题目描述
复合梯形公式是数值积分中最基本、最常用的方法之一。它的核心思想是将积分区间等分成多个小区间,在每个小区间上用梯形面积来近似积分。然而,其收敛速度较慢,截断误差与步长的平方成正比。本题目将详细讲解如何利用 Richardson 外推技术对复合梯形公式进行加速,使其在相同计算量下获得更高的精度。
解题过程循序渐进讲解
第一步:回顾复合梯形公式及其误差
给定一个在区间 [a, b] 上可积的函数 f(x),我们希望计算定积分:
\[I = \int_{a}^{b} f(x) \, dx. \]
将区间 [a, b] 等分为 n 个子区间,步长 h = (b - a) / n,节点为 \(x_i = a + i h, i = 0, 1, \dots, n\)。则复合梯形公式为:
\[T(h) = h \left[ \frac{1}{2} f(a) + f(x_1) + f(x_2) + \dots + f(x_{n-1}) + \frac{1}{2} f(b) \right]. \]
误差分析:当 f(x) 足够光滑时,复合梯形公式的截断误差可表示为关于 h 的渐近展开:
\[I - T(h) = c_1 h^2 + c_2 h^4 + c_3 h^6 + \dots \]
其中系数 \(c_k\) 与 f 的高阶导数在端点处的值有关,且展开中只包含 h 的偶次幂。这个展开是 Richardson 外推能够加速的基础。
第二步:理解 Richardson 外推的基本原理
Richardson 外推是一种利用低阶近似值的误差展开式,通过线性组合来消除低阶误差项,从而得到更高阶近似的方法。针对复合梯形公式,我们有:
\[I = T(h) + c_1 h^2 + c_2 h^4 + c_3 h^6 + \dots \]
如果我们用步长 h 和步长 h/2 分别计算梯形公式近似值 T(h) 和 T(h/2),则:
\[I = T(h) + c_1 h^2 + c_2 h^4 + c_3 h^6 + \dots \\ I = T(h/2) + c_1 (h/2)^2 + c_2 (h/2)^4 + c_3 (h/2)^6 + \dots \]
将第二式乘以 4 再减去第一式,可以消去 \(c_1 h^2\) 项:
\[4I - I = 4T(h/2) - T(h) + (4c_2 (h/2)^4 - c_2 h^4) + \dots \]
整理后得到:
\[I = \frac{4T(h/2) - T(h)}{3} - \frac{1}{4} c_2 h^4 + \dots \]
记:
\[S(h) = \frac{4T(h/2) - T(h)}{3}. \]
这个新的近似 S(h) 的误差从原来的 \(O(h^2)\) 提升到了 \(O(h^4)\)。实际上,S(h) 恰好就是复合辛普森公式的近似值(当区间等分时)。
第三步:推广到多级外推(龙贝格积分法的思想)
我们可以重复这一外推过程。对 S(h) 来说,其误差展开为:
\[I - S(h) = d_1 h^4 + d_2 h^6 + d_3 h^8 + \dots \]
注意这里没有 h^2 项,但可能存在 h^6 项(由于对称性,实际展开中只含 h 的偶次幂,但不同外推级的系数会变化)。如果我们再计算 S(h/2) ,类似地组合:
\[R(h) = \frac{16S(h/2) - S(h)}{15}, \]
则 R(h) 的误差变为 \(O(h^6)\)。这个过程可以一直进行下去,得到误差阶越来越高的近似值。这就是龙贝格积分法的核心思想,但本题目聚焦于复合梯形公式的 Richardson 外推技巧,因此我们重点掌握前两到三级的外推。
第四步:计算步骤示例
以积分 \(I = \int_0^1 e^{-x^2} dx\) 为例,演示如何应用。
- 选择初始步长:取初始步长 \(h_1 = 1\) (即整个区间作为一个梯形),计算:
\[ T(h_1) = 1 \times \left[ \frac{1}{2}f(0) + \frac{1}{2}f(1) \right] = \frac{1}{2}(1 + e^{-1}) \approx 0.6839397206. \]
- 将步长减半:取 \(h_2 = 1/2\),需要计算新增节点 \(x=0.5\) 处的函数值。梯形公式为:
\[ T(h_2) = 0.5 \times \left[ \frac{1}{2}f(0) + f(0.5) + \frac{1}{2}f(1) \right] \\ = 0.5 \times [0.5 \times 1 + e^{-0.25} + 0.5 \times e^{-1}] \approx 0.7313702518. \]
- 进行一次 Richardson 外推(即得到复合辛普森近似):
\[ S(h_2) = \frac{4T(h_2) - T(h_1)}{3} = \frac{4 \times 0.7313702518 - 0.6839397206}{3} \approx 0.7471804289. \]
此时误差阶为 \(O(h^4)\)。
- 再次将步长减半(可选,演示二级外推):取 \(h_3 = 1/4\),计算 \(T(h_3)\) 需要新增节点 0.25, 0.75 处的函数值。
\[ T(h_3) = 0.25 \times [0.5f(0) + f(0.25) + f(0.5) + f(0.75) + 0.5f(1)] \\ = 0.25 \times [0.5 + e^{-0.0625} + e^{-0.25} + e^{-0.5625} + 0.5e^{-1}] \approx 0.7429840978. \]
- 利用 T(h_2) 和 T(h_3) 进行第一级外推得到 S(h_3):
\[ S(h_3) = \frac{4T(h_3) - T(h_2)}{3} \approx \frac{4 \times 0.7429840978 - 0.7313702518}{3} \approx 0.7468553798. \]
- 进行第二级 Richardson 外推,利用 S(h_2) 和 S(h_3):
\[ R(h_3) = \frac{16S(h_3) - S(h_2)}{15} \approx \frac{16 \times 0.7468553798 - 0.7471804289}{15} \approx 0.7468337093. \]
这个 R(h_3) 的误差阶为 \(O(h^6)\)。与精确值 0.7468241329 相比,误差约为 9.6e-6,而原始的 T(h_3) 误差约为 3.8e-3,可见外推显著提高了精度。
第五步:误差与收敛性讨论
Richardson 外推技巧的有效性依赖于误差展开式中存在 h 的幂次项。对于复合梯形公式,由于被积函数足够光滑时误差展开只有 h 的偶次幂,所以外推可以逐步消除低阶项。但在某些情况下(如被积函数具有奇异性或振荡剧烈),误差展开可能不符合形式,外推效果会变差甚至发散。因此,在实际应用中,通常结合自适应策略,在误差展开有效的区间内进行外推。
第六步:算法实现要点
- 从一个较粗的步长开始,逐步二分网格,每次只需计算新增节点的函数值,可以利用之前的结果高效计算新的梯形公式近似。
- 建立一个外推表(类似龙贝格表),第一列为不同步长的梯形近似 T,后续每一列通过 Richardson 外推公式从前一列计算得到,直至满足精度要求或外推收敛。
- 终止条件可以设置为相邻两次外推结果的绝对差或相对差小于给定容差。
通过以上步骤,你可以清晰地看到,Richardson 外推技巧能够在不显著增加计算量的情况下(每次二分只增加约一倍的函数求值,但精度提升多阶),大幅提高复合梯形公式的收敛速度。这是一种简单而强大的数值积分加速技术。