数值积分中的外推法
外推法是一种强大的技术,它能够显著提高数值结果的精度,而无需显著增加计算成本。其核心思想是:通过对一系列粗糙近似的结果进行巧妙的线性组合,来“外推”出更精确的近似值。这种方法在数值积分、微分、微分方程求解等许多计算数学领域都有广泛应用。下面,我们循序渐进地来理解它。
第一步:从误差展开式开始
理解外推法的关键在于认识到许多数值方法的近似误差常常具有一个已知的渐进展开式。以数值积分为例,考虑用复合梯形法则(Trapezoidal Rule)计算定积分 \(I = \int_a^b f(x) dx\)。
假设我们将积分区间 \([a, b]\) 等分为 \(n\) 个子区间,步长为 \(h = (b-a)/n\),得到的近似积分值记为 \(T(h)\)。对于充分光滑的被积函数 \(f(x)\),其误差可以展开为 \(h\) 的偶数次幂级数(欧拉-麦克劳林公式):
\[I = T(h) + c_2 h^2 + c_4 h^4 + c_6 h^6 + \dots \]
其中 \(c_2, c_4, \dots\) 是与 \(f\) 在端点的导数值有关的系数,与 \(h\) 无关。
这个式子告诉我们,\(T(h)\) 的误差主要是 \(O(h^2)\) 量级的。如果我们用不同的步长计算,比如用 \(h\) 和 \(h/2\),就会得到两个不同精度的近似值 \(T(h)\) 和 \(T(h/2)\)。
第二步:构建理查德森外推
理查德森外推是外推法最经典和基本的形式。我们利用上一步的误差展开式。
- 写出两个近似:
- 用步长 \(h\): \(I = T(h) + c_2 h^2 + c_4 h^4 + \dots\)
- 用步长 \(h/2\): \(I = T(h/2) + c_2 (h/2)^2 + c_4 (h/2)^4 + \dots = T(h/2) + \frac{c_2}{4}h^2 + \frac{c_4}{16}h^4 + \dots\)
- 消除主误差项:观察发现,这两个等式中的主误差项(即 \(h^2\) 项)的系数相差4倍。我们可以将第二个等式乘以4,减去第一个等式:
\[ 4I - I = [4T(h/2) - T(h)] + [c_2h^2 - c_2h^2] + [\frac{c_4}{4}h^4 - c_4 h^4] + \dots \]
整理后得到:
\[ I = \frac{4T(h/2) - T(h)}{4 - 1} - \frac{c_4}{4}h^4 + \dots \]
- 得到高阶近似:我们定义一个新的近似:
\[ S(h/2) = \frac{4T(h/2) - T(h)}{3} \]
这个新近似 \(S(h/2)\) 的误差项从原来的 \(O(h^2)\) 变成了 \(O(h^4)\)!我们仅用 \(T(h)\) 和 \(T(h/2)\) 这两个精度较低的结果,通过一个线性组合,就得到了一个精度高得多的结果。这个过程就是一次理查德森外推。在这里,\(S(h/2)\) 实际上等价于复合辛普森法则。
第三步:龙贝格积分——重复外推的过程
龙贝格积分是将理查德森外推法系统化、自动化地应用于梯形法则,以高效计算高精度积分值的算法。
- 生成初始梯形序列:首先,计算一系列步长不断减半的复合梯形近似值。记 \(R_{k,0}\) 为第 \(k\) 次二分后的梯形近似值(\(k=0,1,2,\dots\))。
- \(R_{0,0}\): 区间 \([a,b]\) 不分割的梯形公式(步长 \(h_0 = b-a\))。
- \(R_{1,0}\): 二分一次(步长 \(h_1 = h_0/2\))的梯形公式。
- \(R_{2,0}\): 再二分一次(步长 \(h_2 = h_1/2\))的梯形公式。
- ……
- 建立外推表:我们利用误差展开式 \(I = R_{k,0} + c_2 h_k^2 + c_4 h_k^4 + \dots\),其中 \(h_k = (b-a)/2^k\)。对相邻的两行应用理查德森外推,可以消除 \(h^2\) 项:
\[ R_{k,1} = \frac{4^1 R_{k,0} - R_{k-1,0}}{4^1 - 1} \]
这样得到的序列 \(R_{1,1}, R_{2,1}, R_{3,1}, \dots\) 的误差是 \(O(h^4)\) 量级,它们等价于辛普森法则序列。
- 继续外推:新序列 \(R_{k,1}\) 的误差展开式以 \(h^4\) 为主项。我们可以继续对它们进行外推,以消除 \(h^4\) 项:
\[ R_{k,2} = \frac{4^2 R_{k,1} - R_{k-1,1}}{4^2 - 1} \]
这个新序列 \(R_{k,2}\) 的误差是 \(O(h^6)\) 量级,等价于布尔法则。
- 一般公式与三角表:这个过程可以持续进行。一般的外推公式为:
\[ R_{k,j} = \frac{4^j R_{k, j-1} - R_{k-1, j-1}}{4^j - 1}, \quad k = 1,2,\dots; \quad j = 1,2,\dots, k \]
将所有结果排列成一个下三角阵(龙贝格表):
```
R_{0,0}
R_{1,0} R_{1,1}
R_{2,0} R_{2,1} R_{2,2}
R_{3,0} R_{3,1} R_{3,2} R_{3,3}
... ... ... ... ...
```
表中对角线上的元素 \(R_{0,0}, R_{1,1}, R_{2,2}, \dots\) 精度依次提高。算法通常在这些值的变化小于某个预设容差时停止。
第四步:扩展与应用要点
-
更一般的外推公式:理查德森外推不仅限于 \(h^2\) 展开。如果已知某种近似 \(A(h)\) 满足误差展开式 \(I = A(h) + c_1 h^{p} + c_2 h^{q} + \dots\)(其中 \(p < q < \dots\)),就可以用类似的思想组合 \(A(h)\) 和 \(A(h/t)\)(\(t\) 是步长收缩因子)来消除 \(h^p\) 项。
-
计算优势:外推法的美妙之处在于,它用相对便宜的低阶方法计算结果,通过“后期处理”(线性组合)获得了高阶方法的精度,而无需直接实现复杂的高阶公式。在龙贝格积分中,每次二分梯形序列只需计算新增节点的函数值,计算量增加不多,但通过外推,精度呈指数级(2的幂次)提高。
-
注意事项:
- 光滑性前提:误差的渐进展开式成立要求被积函数(或被近似的对象)足够光滑。如果函数有奇点、间断点或不光滑,展开式可能不成立,外推会失效甚至给出错误结果。
- 舍入误差累积:外推过程涉及数值相减,当外推次数过多时,可能会放大舍入误差,导致精度无法继续提高,甚至下降。因此,外推次数是有限制的。
总结来说,数值积分中的外推法(以理查德森外推和龙贝格积分为代表)是一种基于误差渐进展开的加速收敛技术。它从一个低阶、经济的数值方法序列出发,通过巧妙的线性组合,逐级消去误差展开式中的低阶项,从而高效、系统地生成一系列精度迅速提高的近似值,是用“智慧”换取计算效率的典范。