复化求积方法对提高精度是行之有效的,但在使用求积公式之前必须先给出步长。步长取得太大精度难以保证,步长太小则又会导致计算量的增加,而事先给出一个合适的步长往往是困难的。
实际计算过程中常采用变步长的计算方案,即在步长逐次减半(即步长二分)的过程中,反复利用复化求积公式进行计算,直到所求得的积分值满足精度要求为止。
设将求积区间 [a, b] 分为 n 等分,则一共有 n + 1 个等分点 x k = a + k h x_{k}=a+kh xk=a+kh, h = b − a n h=\frac{b-a}{n} h=nb−a, k = 0 , 1 , . . . , n k=0, 1 ,..., n k=0,1,...,n。这里用 T n T_{n} Tn 表示用复化梯形法求得的积分值,其下表 n 表示等分数。
先考察一个字段 [ x k , x k + 1 ] [x_{k}, x_{k+1}] [xk,xk+1],其中点 x k + 1 2 = 1 2 ( x k + x k + 1 ) x_{k+\frac{1}{2}}=\frac{1}{2}(x_{k}+x_{k+1}) xk+21=21(xk+xk+1),在该子段上二分前后的两个积分值
T
1
=
h
2
[
f
(
x
k
)
+
f
(
x
k
+
1
)
]
T_{1}=\frac{h}{2}[f(x_{k})+f(x_{k+1})]
T1=2h[f(xk)+f(xk+1)]
T
2
=
h
4
[
f
(
x
k
)
+
2
f
(
x
k
+
1
2
)
+
f
(
x
k
+
1
)
]
T_{2}=\frac{h}{4}[f(x_{k})+2f(x_{k+\frac{1}{2}})+f(x_{k+1})]
T2=4h[f(xk)+2f(xk+21)+f(xk+1)]
显然有下列关系
T 2 = 1 2 T 1 + h 2 f ( x k + 1 2 ) T_{2}=\frac{1}{2}T_{1}+\frac{h}{2}f(x_{k+\frac{1}{2}}) T2=21T1+2hf(xk+21)
将这一关系式关于 k 从 0 到 n - 1 累加求和,即可导出下列递推公式
T 2 n = 1 2 T n + h 2 ∑ k = 0 n − 1 f ( x k + 1 2 ) (1) T_{2n}=\frac{1}{2}T_{n}+\frac{h}{2}\sum_{k=0}^{n-1} f(x_{k+\frac{1}{2}})\tag1 T2n=21Tn+2hk=0∑n−1f(xk+21)(1)
需要强调指出的是,上式中的 h = b − a n h=\frac{b-a}{n} h=nb−a 代表二分前的步长,而
x k + 1 2 = a + ( k + 1 2 ) h x_{k+\frac{1}{2}}=a+(k+\frac{1}{2})h xk+21=a+(k+21)h
图 2-2 是变步长梯形法的算法框图。其中 T 1 T_{1} T1 和 T 2 T_{2} T2 分别代表二分前后的积分值。各框的含义是:
- [框 1] 准备初值
- [框 2] 按式(1)求二分后的梯形值
从第一个分点 x = a + h 2 x=a+\frac{h}{2} x=a+2h 出发,取 h 为步长逐步向右跨,即可依次确定式(1)中的各个分点。图 2-2 将所得到的分点暂存于单元 x 中。
- [框 3] 控制精度
- [框 4] 修改步长
最终取二分后的积分值 T 2 T_{2} T2 作为所求的结果。
例:用变步长方法计算
I = ∫ 0 1 sin x x d x I=\int_{0}^{1}\frac{\sin x}{x}dx I=∫01xsinxdx
解:先对整个区间 [0, 1] 用梯形公式。 f ( x ) = sin x x f(x)=\frac{\sin x}{x} f(x)=xsinx 的数据参见下表。
x | f(x) |
---|---|
0 | 1.000 000 0 |
1/8 | 0.997 397 8 |
1/4 | 0.989 615 8 |
3/8 | 0.916 726 7 |
1/2 | 0.958 851 0 |
5/8 | 0.936 155 6 |
3/4 | 0.908 851 6 |
7/8 | 0.877 192 5 |
1 | 0.841 470 9 |
由于 f ( 0 ) = 1 f(0)=1 f(0)=1, f ( 1 ) = 0.8414709 f(1)=0.8414709 f(1)=0.8414709,故
T 1 = 1 2 [ f ( 0 ) + f ( 1 ) ] = 0.9207355 T_{1}=\frac{1}{2}[f(0)+f(1)]=0.9207355 T1=21[f(0)+f(1)]=0.9207355
然后将区间二分,由于 f ( 1 2 ) = 0.9588510 f(\frac{1}{2})=0.9588510 f(21)=0.9588510,利用递推公式(1)得
T 2 = 1 2 T 1 + 1 2 f ( 1 2 ) = 0.9397933 T_{2}=\frac{1}{2}T_{1}+\frac{1}{2}f(\frac{1}{2})=0.9397933 T2=21T1+21f(21)=0.9397933
再二分一次,并计算新分点上的函数值
f ( 1 4 ) = 0.9896158 f(\frac{1}{4})=0.9896158 f(41)=0.9896158
f ( 3 4 ) = 0.9088516 f(\frac{3}{4})=0.9088516 f(43)=0.9088516
再用式(1)求得
T 4 = 1 2 T 2 + 1 4 [ f ( 1 4 ) + f ( 3 4 ) ] = 0.9445135 T_{4}=\frac{1}{2}T_{2}+\frac{1}{4}[f(\frac{1}{4})+f(\frac{3}{4})]=0.9445135 T4=21T2+41[f(41)+f(43)]=0.9445135
这样不断二分下去,计算结果如下表所示(表中 k 代表二分次数,区间等分数 n = 2k)。积分的准确值为 0.9460831,用变步长方法二分 10 次得到了这个结果。
k | Tn |
---|---|
0 | 0.920 735 5 |
1 | 0.939 793 3 |
2 | 0.944 513 5 |
3 | 0.945 690 9 |
4 | 0.945 985 0 |
5 | 0.946 058 6 |
6 | 0.946 076 9 |
7 | 0.946 081 5 |
8 | 0.946 082 7 |
9 | 0.946 083 0 |
10 | 0.946 083 1 |
运行示例:
程序源码:
#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;
int count = 0;
double T2;
double h = b - a;
double T1 = h / 2 * (1 + f(b));
do
{
printf("第 %d 次二分步长,二分后的积分值 T2 = %.7f\n", ++count, T1);
// sum 为各分点的函数值的和
double sum = 0;
// 分点的值
double x = a + h / 2;
while (x < b)
{
// 将分点往后区间范围内的分点函数值求和
sum += f(x);
// 分点往后移动(在区间上限范围内)
x += h;
}
// 求得步长二分后的积分值
T2 = T1 / 2 + h / 2 * sum;
h /= 2;
double temp = T1;
T1 = T2;
T2 = temp;
} while (abs(T2 - T1) - accuracy > 0);
return 0;
}