最后更新: 2026-06-27 17:05 查看: 3 次纠错评论(0) 课件开VIP

外推法

7．4．1 理查森外推
计算数学中最重要的思想之一是我们可以从一些近似中获取信息，使用它们来估计近似误差并生成更好的改进方法。在本节中，我们将对其中一些递推方法进行介绍。

考虑一个通用数值积分关系 $I_n(f)$ 来近似给定一个积分 $I(f)$ 。假定该积分的误差关系为

$$
\begin{equation*}
I(f)-I_n(f) \approx C h^{-m} . \tag{7.20}
\end{equation*}
$$

比如，在梯形公式中 $m=2$ ，中点公式中 $m=2$ 和辛普森公式中 $m=4$ 都是正确的。基于前面的假设，有

$$
I-I_{2 n} \approx C(2 h)^{-m}=2^{-m}\left(C h^{-m}\right) \approx 2^{-m}\left(I-I_n\right)
$$

然后我们可以求解近似等式

$$
I \approx \frac{2^m I_{2 n}-I_n}{2^m-1}
$$

从而我们得到如下的理查森（Richardson）外推

$$
R_{2 n}=\frac{2^m I_{2 n}-I_n}{2^m-1}
$$

并且有可计算误差估计

$$
\begin{equation*}
E_{2 n}=R_{2 n}-I_{2 n}=\frac{I_{2 n}-I_n}{2^m-1} \tag{7.21}
\end{equation*}
$$

这样一来，将外推值 $R_{2 n}$ 作为精确值的近似比直接的 $I_{2 n}$ 更好。比如利用表7．7中的数据我们得到

![图片](/uploads/2026-06/0afadf.jpg)
 
$$
S_8(f)=0.746826120527, \quad S_{16}(f)=0.746824257436 .
$$

于是

$$
R_{16}(f)=0.7468241335, \quad E_{16}(f)=0.12421 \times 10^{-6} .
$$

下面是利用理查森外推法计算 $\int_0^1 \mathrm{e}^x \mathrm{~d} x$ 的例子。表7．7显示了基于（7．20）和（7．21）的外推近似值和误差估计以及实际误差．随着网格尺寸减半，误差 $I(f)-R_n(f)$ 减少了大约 16 倍，这表明梯形公式的误差已通过外推法从 $O\left(h^2\right)$ 改进为 $O\left(h^4\right)$ ．

例 7.8 考虑积分

$$
I(f)=\int_0^1 \mathrm{e}^{-x^3} \mathrm{~d} x
$$

对于该积分函数，我们无法得到精确值，这样在数值计算过程中我们无法将近似值与精确解进行比较．但是，如果我们用辛普森公式计算该积分值，则可以得出

$$
S_{16}(f)=0.807512351889, \quad S_{32}(f)=0.807511254956 .
$$

误差估计为

$$
E_{32}=0.731 \times 10^{-7} .
$$

由于 $E_{32}<10^{-7}$ 很小，因此我们有理由认为实际计算的误差也很小．其实利用外推值

$$
R_{32}=0.807511182
$$

会更加准确．所以理查森外推可以用来估计计算误差．
7．4．2 龙贝格积分法
首先不加证明地直接引用如下的欧拉－麦克劳林（Euler－Maclaurin）公式．

定理7．7（欧拉－麦克劳林公式）如果 $f$ 具有足够的可微性，则对于任何 $N>0$ ，都有一组常数 $c_k, \quad 1 \leqslant k \leqslant N+1$ ，使得对于 $\xi \in[a, b]$ ，梯形公式中的误差满足

$$
\begin{equation*}
I(f)-T_n(f)=\gamma_1 h^2+\gamma_2 h^4+\cdots+\gamma_N h^{2 N}+c_{N+1}(b-a) h^{2 N+2} f^{(2 N+2)}(\xi), \tag{7.22}
\end{equation*}
$$

其中 $\gamma_k=c_k\left(f^{(2 k-1)}(b)-f^{(2 k-1)}(a)\right)$ ．
欧拉—麦克劳林公式的意义在于，它使我们可以将梯形公式中的误差写为网格间距 $h$ 的幂次数，然后使用该序列来推导越来越精确的新积分关系。注意，只要 $f(x)$ 可微，就可以进行误差的级数展开。欧拉－麦克劳林公式的另一个结果是，当将梯形公式应用到周期函数时显示出来的精度特别高。

下面把误差扩展，以构建更精确的积分关系，只需将（7．22）子区间个数扩大两倍，即将 $h$ 替换为 $h / 2$ ，

$$
\begin{equation*}
I(f)-T_{2 n}(f)=\frac{1}{4} \gamma_1 h^2+\frac{1}{16} \gamma_2 h^4+\cdots+\frac{1}{4^N} \gamma_N h^{2 N}+O\left(h^{2 N+2}\right) . \tag{7.23}
\end{equation*}
$$

如果将（7．23）乘以 4 ，再减去（7．22），然后求出 $I(f)$ ，我们得到

$$
\begin{equation*}
I(f)-\frac{4 T_{2 n}(f)-T_n(f)}{3}=b_2 h^4+b_3 h^6+\cdots+b_N h^{2 N}+O\left(h^{2 N+2}\right), \tag{7.24}
\end{equation*}
$$

其中 $b_k=-\frac{1}{3}\left(1-4^{1-k}\right) c_k\left(f^{(2 k-1)}(b)-f^{(2 k-1)}(a)\right)$ ．

注意，这里值 $\frac{1}{3}\left(4 T_{2 n}(f)-T_n(f)\right)$ 恰好是理查森对于 $T_n(f)$ 的外推值（7．21）．表达式（7．24）表明它是 $O\left(h^4\right)$ 的．

我们将（7．24）改写为

$$
\begin{equation*}
I(f)-\frac{4 T_{2 n}(f)-T_n(f)}{3}=b_2 h^4+b_3 h^6+\cdots+b_N h^{2 N}+O\left(h^{2 N+2}\right) . \tag{7.25}
\end{equation*}
$$

为此，我们还有

$$
\begin{equation*}
I(f)-R_{4 n}(f)=\frac{1}{16} b_2 h^4+\frac{1}{64} b_3 h^6+\cdots+\frac{1}{4^N} b_N h^{2 N}+O\left(h^{2 N+2}\right) . \tag{7.26}
\end{equation*}
$$

将（7．26）乘以 16，再减去（7．25），得到

$$
I(f)-\frac{16 R_{4 n}-R_{2 n}}{15}=a_3 h^6+\cdots+a_k h^{2 k}+\cdots+a_N h^{2 N}+O\left(h^{2 N+2}\right),
$$

其中 $a_k=-\frac{1}{15}\left(1-4^{2-k}\right) b_k$ ．因此，$\frac{1}{15}\left(16 R_{4 n}-R_{2 n}\right)$ 可以看作是积分的另一种外推近似．
外推过程的每个步骤都会产生一个新的积分组合方法，该方法比前一个方法更准确．可以将该过程进行系统化，以产生称为龙贝格（Romberg）积分的算法．

为更进一步讨论，我们需要定义一些符号．令 $T_n^{(0)}(f)$ 表示梯形公式值，使其成为三角形数组的第一列。我们用 $T_n^{(1)}(f)$ 表示第二列。根据理查森外推公式从梯形值中计算这些值

$$
T_{2 n}^{(1)}(f)=R_{2 n}(f)=\frac{4 T_{2 n}^{(0)}(f)-T_n^{(0)}(f)}{3} .
$$

请注意，下标意味着第二列将比第一列少一个数据．通常，根据公式从上一列中计算出每一列

$$
\begin{equation*}
T_{2 n}^{(j+1)}(f)=\frac{4^{j+1} T_{2 n}^{(j)}(f)-T_n^{(j)}(f)}{4^{j+1}-1} . \tag{7.27}
\end{equation*}
$$

这样一来，计算序列可以列为如下形式：

$$
\begin{array}{cccccc}
T_n^{(0)}(f) & & & & \\
T_{2 n}^{(0)}(f) & T_{2 n}^{(1)}(f) & & & \\
T_{4 n}^{(0)}(f) & T_{4 n}^{(1)}(f) & T_{4 n}^{(2)}(f) & & \\
\vdots & \vdots & \vdots & \ddots & \\
T_{2^k n}^{(0)}(f) & T_{2^k n}^{(1)}(f) & T_{2^k n}^{(2)}(f) & \ldots & T_{2^k n}^{(k)}(f)
\end{array}
$$
如果被积数足够光滑，则数组行中的每一步都将消除误差扩展的 $h^2$ 幂次，即第一列的误差为 $O\left(h^2\right)$ ，第二列误差为 $O\left(h^4\right)$ ，第三列误差为 $O\left(h^6\right)$ ，以此类推。同时，沿着阵列向下移动会使每一行的 $h$ 减少 2 倍．如果我们从 $h=1 / 2$ 开始，那么在第六行中，有 $h=1 / 64$ ．龙贝格积分选择阵列的对角元素来近似，这种近似精度非常高．一般有如下结果。

定理7．8（龙贝格积分）假设 $f(x)$ 在 $[a, b]$ 上可微，令 $\theta_k(f)$ 为龙贝格数组的 $k$ 个对角元素 $\theta_k(f)=T_{2^k n}^{(k)}(f)$ 。如果初始梯形公式 $T_n^{(0)}(f)$ 的网格尺寸为 $h$ ，则

$$
I(f)-\theta_k(f)=O\left(4^{-k} h^{2 k+2}\right)
$$

另外，在龙贝格积分产生准确近似值的同时，它也可以用于估计误差．假设被积函数足够光滑，因为对于龙贝格数组中的每个计算项，我们都可以使用理查森外推法来计算误差的估计值

$$
I(f)-T_{2 n}^{(k)}=E_{2 n}^{(k)}=\frac{T_{2 n}^{(k)}-T_n^{(k)}}{4^k-1}
$$

所以当估计的误差足够小时，便可以停止龙贝格过程．
值得说明的是，龙贝格中的大部分计算工作都涉及计算第一列（梯形公式值），因为由（7．27）完成的后续列的计算仅涉及一些简单的操作。因此，龙贝格的计算效率取决于我们快速递归计算梯形公式值的能力，即比较简单地从 $T_n^{(0)}(f)$ 计算 $T_{2 n}^{(0)}(f)$ 。一般的计算方法如下定理所示。

定理7．9 令 $T_n(f)$ 表示在给定区间 $[a, b]$ 内，通过使用 $n$ 个子区间的梯形公式计算结果，均匀网格间距 $h=(b-a) / n$ ．然后，我们可以使用如下公式来计算两倍的子区间上的梯形公式 $T_{2 n}(f)$ ，

$$
\begin{equation*}
T_{2 n}(f)=\frac{1}{2} T_n(f)+\left(\frac{b-a}{2 n}\right) \sum_{j=1}^n f(a+(2 j-1)(b-a) /(2 n)) . \tag{7.28}
\end{equation*}
$$

证明 关键步骤是认识到我们可以为 $n$ 个子区间编写梯形公式

$$
T_n(f)=\left(\frac{b-a}{n}\right)\left(\frac{1}{2} f(a)+\frac{1}{2} f(b)+\sum_{i=1}^{n-1} f(a+i(b-a) / n)\right) .
$$

因此， $2 n$ 个子区间的公式是

$$
T_{2 n}(f)=\left(\frac{b-a}{2 n}\right)\left(\frac{1}{2} f(a)+\frac{1}{2} f(b)+\sum_{i=1}^{2 n-1} f(a+i(b-a) /(2 n))\right)
$$

现在计算

$$
T_{2 n}(f)=\left(\frac{b-a}{2 n}\right)\left(\frac{1}{2} f(a)+\frac{1}{2} f(b)+\sum_{i=1}^{2 n-1} f(a+i(b-a) /(2 n))\right)
$$

$$
\begin{aligned}
& =\frac{1}{2}\left(\frac{b-a}{n}\right)\left(\frac{1}{2} f(a)+\frac{1}{2} f(b)+\sum_{i=2,4,6, \cdots}^{2 n-1} f(a+i(b-a) /(2 n))+\sum_{i=1,3,5}^{2 n-1} f(a+i(b-a) /(2 n))\right) \\
& =\frac{1}{2}\left(\frac{b-a}{n}\right)\left(\frac{1}{2} f(a)+\frac{1}{2} f(b)+\sum_{j=1}^{n-1} f(a+2 j(b-a) /(2 n))+\sum_{j=1}^n f(a+(2 j-1)(b-a) /(2 n))\right) \\
& =\frac{1}{2}\left(\frac{b-a}{n}\right)\left(\left[\frac{1}{2} f(a)+\frac{1}{2} f(b)+\sum_{j=1}^{n-1} f(a+j(b-a) / n)\right]+\sum_{j=1}^n f(a+(2 j-1)(b-a) /(2 n))\right) \\
& =\frac{1}{2} T_n(f)+\left(\frac{b-a}{2 n}\right)\left(\sum_{j=1}^n f(a+(2 j-1)(b-a) /(2 n))\right) .
\end{aligned}
$$

这样就完成了证明．
这个结果的重点是我们可以使用最少数量的函数求值来计算龙贝格数组的整个第一列。 $T_n^{(0)}(f)$ 是从头开始计算的，将要求我们重新计算 $f(x)$ 的 $n$ 个数值。

这个定理告诉我们，可以使用不超过普通梯形公式的函数求值来计算 $T_n^{(0)}(f)$ ，但是我们可以获得更高的准确性。表7．8显示了龙贝格积分在我们的标准示例中的应用，该示例在区间 $[0,1]$ 上对 $\mathrm{e}^x$ 进行积分，这里的误差迅速衰减．

![图片](/uploads/2026-06/fcab95.jpg)
 
 算法7．1 龙贝格积分
1．计算 $T_n^{(0)}(f)=T_n(f)$ ，这是初始梯形公式计算，使用 $n$ 个子区间（通常 $n=1$ ，但不是必需的）．

2．对于从 1 到 $N$ 的 $k$ ，
（1）计算 $T_{2^k n}^{(0)}(f)=T_{2^k n}(f)$ ，这是龙贝格数组新行上的第一个数值．
（2）在行上外推：对于 $j$ ，从 0 到 $k-1$ 计算

$$
T_{2^k n}^{(j+1)}(f)=\left(4^{j+1} T_{2^k n}^{(j)}(f)-T_{2^{k-1} n}^{(j)}(f)\right) /\left(4^{j+1}-1\right) .
$$

实际编程计算中我们不需要存储整个数组，只需存储数组的一行，而在计算新值时将其覆盖即可。

例 7.9 利用外推法求 $f(x)=\left(1+x^4\right)^{-1},[a, b]=[0,1]$ 的积分．为了以最少的工作量，给出积分的基本梯形公式值（表7．9），即 $T_n^{(0)}(f)=T_n(f)$ 近似

$$
I(f)=\int_0^1 \frac{\mathrm{~d} x}{1+x^4}=0.8669729871 .
$$

![图片](/uploads/2026-06/a86eb1.jpg)
 
 第一次外推产生第一个龙贝格值

$$
T_2^{(1)}(f)=\frac{4 T_2^{(0)}(f)-T_1^{(0)}(f)}{3}=0.8774509804
$$

需要两个推断才能得出第二个龙贝格值

$$
T_4^{(1)}(f)=\frac{4 T_4^{(0)}(f)-T_2^{(0)}(f)}{3}=0.8671137006,
$$

$$
T_4^{(2)}(f)=\frac{16 T_4^{(1)}(f)-T_2^{(1)}(f)}{15}=0.86642454864
$$

$T_4^{(2)}(f)$ 误差估计由下式给出

$$
E_4^{(1)}(\theta)=\left|I(f)-T_4^{(2)}(f)\right|=0.0005484385 .
$$

因此，我们确信龙贝格值 $T_4^{(2)}(f)$ 可以精确到 $10^{-3}$ 以内．

科数题库(单机版) 会议室预约系统(book)

今日还可看 0 篇未注册用户每天查看4篇，注册用户每天8篇，开通VIP会员无限制查看。

免费注册

《高等数学》难点解析

高数教程 泰勒公式切线与法线切平面与法平面驻点·拐点·极值点·零点间断点渐进线瑕积分欧拉方程伯努利方程 Abel 收敛定理偏导数的几何意义偏导数的几何意义梯度数量场与向量场多元函数极值拉格朗日算子通量与散度环流量与旋度格林公式高斯公式斯托克斯公式三大公式比较傅里叶级数极坐标微元点法式方程变上限定积分 X型计算面积 Y型计算面积微分的意义渐近线间断点 y''+py'+qy=f(x)方程高斯黎曼傅里叶变换(复数) 拉普拉斯变换(复数)

《线性代数》难点解析

线代教程 近世代数对数学的整体思考线性的意义矩阵乘法(列视角) 矩阵乘法(行视角) 矩阵左乘矩阵右乘逆矩阵求解方程组阶梯形矩阵的求法方程组解的判定四阶行列式的计算线性变换的意义线性空间向量组的等价线性空间的几何意义基础解系的求法施密特正交化特征值与特征向量的意义矩阵相似的几何意义矩阵可对角化的理解秩的意义(向量版) 秩的意义(方程版) 二次型的意义

《概率论与数理统计》难点解析

概率教程 置信区间与上a分位数概率中的“取”与“放” 贝叶斯公式全概率公式泊松分布指数分布伽玛分布二维密度图的意义卷积的意义相关系数的意义 k阶矩是与矩母函数卡方分布的作用单正态区间估计理解假设检验理解切比雪夫不等式中心极限定理

上一篇：高斯积分

下一篇：自适应积分

本文对您是否有用？有用 (0) 无用 (0)