梯形法 的算法简单,但精度低,收敛速度缓慢。如何提高收敛速度以节省计算量,自然是人们极为关心的课题。
根据梯形法的误差公式
I − T n ≈ − h 2 12 [ f ′ ( b ) − f ′ ( a ) ] I-T_{n}\approx -\frac{h^{2}}{12}[f'(b)-f'(a)] I−Tn≈−12h2[f′(b)−f′(a)]
积分值 T n T_{n} Tn 的截断误差大致与 h2 成正比,因此步长减半后误差将减至 1 4 \frac{1}{4} 41,即有
I − T 2 n I − T n ≈ 1 4 \frac{I-T_{2n}}{I-T_{n}}\approx \frac{1}{4} I−TnI−T2n≈41
将上式移项整理,知
I − T 2 n ≈ 1 3 ( T 2 n − T n ) (1) I-T_{2n}\approx \frac{1}{3}(T_{2n}-T_{n})\tag1 I−T2n≈31(T2n−Tn)(1)
由此可见,只要二分前后两个积分值 T n T_{n} Tn 与 T 2 n T_{2n} T2n 相当接近,就可以保证计算结果 T 2 n T_{2n} T2n 的误差很小。
按式(1),积分值 T 2 n T_{2n} T2n 的误差大致等于 1 3 ( T 2 n − T n ) \frac{1}{3}(T_{2n}-T{n}) 31(T2n−Tn),如果用这个误差值作为 T 2 n T_{2n} T2n 的一种补偿,可以期望所得的
T ‾ = T 2 n + 1 3 ( T 2 n − T n ) = 4 3 T 2 n − 1 3 T n (2) \overline T=T_{2n}+\frac{1}{3}(T_{2n}-T_{n})=\frac{4}{3}T_{2n}-\frac{1}{3}T_{n}\tag2 T=T2n+31(T2n−Tn)=34T2n−31Tn(2)
应当是更好的结果。
再考察 此例,所求得的两个梯形值 T 4 = 0.9445135 T_{4}=0.944 513 5 T4=0.9445135 和 T 8 = 0.9456909 T_{8}=0.945 690 9 T8=0.9456909 的精度都很差(与准确值 0.9460831 比较,它们只有二三位有效数字),但如果将它们按式(2)做线性组合,则新的近似值
T ‾ = 4 3 T 8 − 1 3 T 4 = 0.9460834 \overline T=\frac{4}{3}T_{8}-\frac{1}{3}T_{4}=0.946 083 4 T=34T8−31T4=0.9460834
却有六位有效数字。
按式(2)组合得到的近似值 T ‾ \overline T T,其实质究竟是什么呢?直接验证易知
S n = 4 3 T 2 n − 1 3 T n (3) S_{n}=\frac{4}{3}T_{2n}-\frac{1}{3}T_{n}\tag3 Sn=34T2n−31Tn(3)
这就是说,用梯形法二分前后两个积分值 T n T_{n} Tn 与 T 2 n T_{2n} T2n 按式(2)作线性组合,结果得到辛普生的积分值 S n S_{n} Sn。
再考察辛普生法。按式
I − S n ≈ − 1 180 ( h 2 ) 4 [ f ′ ′ ′ ( b ) − f ′ ′ ′ ( a ) ] I-S_{n}\approx -\frac{1}{180}(\frac{h}{2})^{4}[f'''(b)-f'''(a)] I−Sn≈−1801(2h)4[f′′′(b)−f′′′(a)]
其截断误差与 h 4 h_{4} h4 成正比。因此,若将步长折半,则误差相应地减少至 1/16,即有
I − S 2 n I − S n ≈ 1 16 \frac{I-S_{2n}}{I-S_{n}}\approx \frac{1}{16} I−SnI−S2n≈161
由此得
I ≈ 16 15 S 2 n − 1 15 S n I\approx \frac{16}{15}S_{2n}-\frac{1}{15}S_{n} I≈1516S2n−151Sn
不难验证,上式右端的值其实等于 C n C_{n} Cn,就是说,用辛普生法二分前后的两个积分值 S n S_{n} Sn 与 S 2 n S_{2n} S2n,按上式再作线性组合,结果得到柯特斯法的积分值 C n C_{n} Cn,即有
C n ≈ 16 15 S 2 n − 1 15 S n (4) C_{n}\approx \frac{16}{15}S_{2n}-\frac{1}{15}S_{n}\tag4 Cn≈1516S2n−151Sn(4)
重复同样的手续,依据柯特斯的误差公式
I − C n ≈ − 2 945 ( h 4 ) 6 [ f ( 5 ) ( b ) − f ( 5 ) ( a ) ] I-C_{n}\approx -\frac{2}{945}(\frac{h}{4})^{6}[f^{(5)}(b)-f^{(5)}(a)] I−Cn≈−9452(4h)6[f(5)(b)−f(5)(a)]
可进一步导出下列龙贝格公式
R n = 64 63 C 2 n − 1 63 C n (5) R_{n}=\frac{64}{63}C_{2n}-\frac{1}{63}C_{n}\tag5 Rn=6364C2n−631Cn(5)
应当注意,龙贝格公式(5)已经不属于牛顿 - 柯特斯公式的范畴。
在步长二分的过程中运行公式(3)~(5)加工三次,就能将粗糙的积分值 T n T_{n} Tn 逐步加工成精度更高的龙贝格值 R n R_{n} Rn,或者说,将收敛缓慢的梯形值序列 T n T_{n} Tn 加工成收敛迅速的龙贝格值序列 R n R_{n} Rn,这种加速算法称为龙贝格算法,其加工流程如图 2-3 所示:
龙贝格算法的程序框图如图 2-4 所示:
例:基于 复化梯形的递推公式的变步长算法求积分 中的例子,用龙贝格算法计算
I = ∫ 0 1 sin x x d x I=\int_{0}^{1}\frac{\sin x}{x}dx I=∫01xsinxdx
得到的梯形值,计算结果如下表所示,表中 k 代表二分次数。
k | T2k | S2k-1 | C2k-2 | R2k-3 |
---|---|---|---|---|
0 | 0.920 135 5 | |||
1 | 0.939 793 3 | 0.946 145 9 | ||
2 | 0.944 513 5 | 0.946 086 9 | 0.946 083 0 | |
3 | 0.945 690 9 | 0.946 083 4 | 0.946 083 1 | 0.946 083 1 |
我们看到,这里用二分 3 次的数据(它们的精度都很低,只有二三位有效数字),通过三次加速获得了例中需要二分 10 次才能得到的结果,而加速过程的计算量(仅仅做几次四则运算而不需要求函数值)可以忽略不计,可见龙贝格加速过程的效果是极其显著的。
运行示例:
程序源码:
#include <iostream>
#include <cmath>
using namespace std;
double f(double x)
{
return sin(x) / x;
}
int main(void)
{
double a, b;
cout << "请输入积分区间:";
cin >> a >> b;
double accuracy;
cout << "请输入精度:";
cin >> accuracy;
// 步长
double h = b - a;
// T1:二分前的梯形法积分值
double T1 = h / 2 * (1 + f(b));
// T2:二分后的梯形法积分值
double T2 = 0;
cout << "T" << pow(2, 0) << " = " << T1 << endl;
// 对 T1 与 T2 加权平均,得辛普森积分值
double S1, S2;
// 对 S1 与 S2 加权平均,得柯特斯积分值
double C1, C2;
// 对 C1 与 C2 加权平均,得龙贝格积分值
double R1, R2;
// flag 作为循环控制标志性变量
int flag = 1;
// 记录二分次数
int k = 1;
while (flag == 1)
{
// 各分点的函数值和
double sum = 0;
// 分点值
double x = a + h / 2;
// 在区间上限范围内求各分点的函数值和
while (x < b)
{
sum += f(x);
x += h;
}
// 计算梯形序列得下一个二分结果
T2 = T1 / 2 + h / 2 * sum;
cout << "T" << pow(2, k) << " = " << T2 << endl;
// 线性组合外推至辛普生
S2 = T2 + 1.0 / 3 * (T2 - T1);
cout << "S" << pow(2, k - 1) << " = " << S2 << endl;
if (k == 1)
{
// 至少外推 2 次得出 S1,S2
k++;
h /= 2;
T1 = T2;
S1 = S2;
continue;
}
else
{
// 线性组合外推值柯特斯
C2 = S2 + 1.0 / 15 * (S2 - S1);
cout << "C" << pow(2, k - 2) << " = " << C2 << endl;
if (k == 2)
{
// 至少外推 3 次得出C1,C2
C1 = C2;
k++;
h /= 2;
T1 = T2;
S1 = S2;
continue;
}
else
{
// 线性组合外推至龙贝格
R2 = C2 + 1.0 / 63 * (C2 - C1);
cout << "R" << pow(2, k - 3) << " = " << R2 << endl;
if (k == 3)
{
// 至少外推 4 次得出 R1, R2
R1 = R2;
C1 = C2;
k++;
h /= 2;
T1 = T2;
S1 = S2;
continue;
}
else if (abs(R2 - R1) >= accuracy)
{
// 精度仍然不符合要求,继续二分步长、继续外推
R1 = R2;
C1 = C2;
k = k + 1;
h = h / 2;
T1 = T2;
S1 = S2;
}
else
{
// 精度符合要求,修改 flag 为 0,跳出 while 循环
flag = 0;
cout << "龙贝格算法求得数值积分结果为:" << R2 << endl;
}
}
}
}
return 0;
}