复合梯形公式的 Richardson 外推加速收敛技巧
题目描述
复合梯形公式是一种基本的数值积分方法,通过将积分区间等分为若干子区间并在每个子区间上应用梯形公式来近似积分值。然而,其收敛速度仅为 \(O(h^2)\),其中 \(h\) 是子区间宽度。为了提高精度,通常需要减小 \(h\),但这会增加计算量。本题目要求利用 Richardson 外推 技术,对复合梯形公式的近似结果进行加速,从而在不显著增加计算成本的前提下获得更高精度的积分值。具体任务是:给定一个函数 \(f(x)\) 在区间 \([a, b]\) 上的积分 \(I = \int_a^b f(x) \, dx\),首先用不同步长(如 \(h, h/2, h/4, \dots\))的复合梯形公式计算一系列近似值,然后通过 Richardson 外推构造一个更高阶的近似序列,最终得到更精确的积分估计。
解题过程
步骤 1:理解复合梯形公式及其误差结构
- 复合梯形公式:将区间 \([a, b]\) 等分为 \(n\) 个子区间,每个子区间宽度 \(h = (b - a)/n\),节点为 \(x_i = a + ih\)(\(i = 0, 1, \dots, n\))。积分近似为:
\[ T(h) = h \left[ \frac{1}{2}f(a) + \sum_{i=1}^{n-1} f(x_i) + \frac{1}{2}f(b) \right]. \]
- 误差分析:复合梯形公式的截断误差可表示为 \(h\) 的偶次幂展开(当 \(f\) 足够光滑时):
\[ I = T(h) + c_2 h^2 + c_4 h^4 + c_6 h^6 + \cdots, \]
其中 \(c_2, c_4, \dots\) 是与 \(f\) 的导数在端点处有关的常数。注意误差中只含 \(h\) 的偶次幂,这是因为梯形公式具有对称性。
步骤 2:Richardson 外推的基本思想
Richardson 外推是一种利用低精度近似序列来消去误差主项、从而得到更高精度近似的技术。
- 假设我们有两个不同步长的近似值 \(T(h)\) 和 \(T(h/2)\)。根据误差展开:
\[ \begin{aligned} I &= T(h) + c_2 h^2 + c_4 h^4 + \cdots, \\ I &= T(h/2) + c_2 (h/2)^2 + c_4 (h/4)^4 + \cdots. \end{aligned} \]
- 将第二个式子乘以 4 并减去第一个式子,可以消去 \(h^2\) 项:
\[ 4I - I = 4T(h/2) - T(h) + \text{高阶误差项}. \]
整理得:
\[ I = \frac{4T(h/2) - T(h)}{3} - \frac{c_4}{4} h^4 + \cdots. \]
定义新的近似:
\[ S(h) = \frac{4T(h/2) - T(h)}{3}. \]
此时 \(S(h)\) 的误差变为 \(O(h^4)\),精度比 \(T(h)\) 提高两阶。
步骤 3:构造外推表格(Richardson 表)
为了系统化外推过程,可以构建一个三角形表格,其中每一行对应不同步长的梯形公式结果,每一列通过外推消去低阶误差项。
- 记 \(R_{k,0} = T\left(\frac{b-a}{2^k}\right)\),即第 \(k\) 行、第 0 列是步长为 \((b-a)/2^k\) 的复合梯形近似。
- 外推递推公式为:
\[ R_{k,m} = \frac{4^m R_{k,m-1} - R_{k-1,m-1}}{4^m - 1}, \quad m = 1, 2, \dots, k. \]
这里 \(m\) 表示外推次数,每外推一次误差阶提高 2(因为消去了一个 \(h^{2m}\) 项)。
3. 表格形式如下(以 \(k=0,1,2\) 为例):
R_{0,0} = T(h)
R_{1,0} = T(h/2) R_{1,1} = (4R_{1,0} - R_{0,0})/3
R_{2,0} = T(h/4) R_{2,1} = (4R_{2,0} - R_{1,0})/3 R_{2,2} = (16R_{2,1} - R_{1,1})/15
其中 \(R_{2,2}\) 的误差为 \(O(h^6)\)。
步骤 4:算法实现步骤
- 输入:被积函数 \(f\),积分区间 \([a, b]\),初始步数 \(n_0\)(如 \(n_0=1\) 对应步长 \(h=b-a\)),外推次数上限 \(M\)。
- 计算第一列(原始梯形公式):
- 对于 \(k=0,1,\dots,M\),计算 \(R_{k,0} = T\left(\frac{b-a}{2^k}\right)\)。注意:可以利用梯形公式的递推关系高效计算,即当步长减半时,只需计算新增节点处的函数值,而不必重新计算所有节点。
- 逐列外推:
对于 \(m=1\) 到 \(M\),对于 \(k=m\) 到 \(M\),计算:
\[ R_{k,m} = R_{k,m-1} + \frac{R_{k,m-1} - R_{k-1,m-1}}{4^m - 1}. \]
(这是递推公式的另一种等价形式,数值上更稳定。)
4. 输出:表格对角线元素 \(R_{k,k}\) 是经过 \(k\) 次外推的结果,误差为 \(O(h^{2k+2})\)。通常取最后一个 \(R_{M,M}\) 作为最终积分近似。
步骤 5:数值稳定性与终止条件
- 当 \(|R_{k,k} - R_{k-1,k-1}| < \varepsilon\)(预设精度)时,可停止外推。
- 由于外推会放大舍入误差,外推次数不宜过多(通常 \(M \leq 10\)),否则高阶结果可能因舍入误差而失真。
- 此方法本质上是 Romberg 积分法 的基础,但 Romberg 积分法通常与梯形公式的递推计算结合,实现更高效。
总结
通过 Richardson 外推,我们可以将复合梯形公式的低阶近似序列转化为高阶近似,显著加速收敛。这种方法不仅适用于梯形公式,也适用于其他具有已知误差展开式的数值方法。在实际应用中,只需逐步加密网格计算梯形值,并利用外推表格不断改进结果,直至满足精度要求。