[计算机数值分析]复化梯形的递推公式的变步长算法求积分

Spring-_-Bear 的 CSDN 博客导航

---

复化求积方法对提高精度是行之有效的，但在使用求积公式之前必须先给出步长。步长取得太大精度难以保证，步长太小则又会导致计算量的增加，而事先给出一个合适的步长往往是困难的。

实际计算过程中常采用变步长的计算方案，即在步长逐次减半（即步长二分）的过程中，反复利用复化求积公式进行计算，直到所求得的积分值满足精度要求为止。

设将求积区间 [a, b] 分为 n 等分，则一共有 n + 1 个等分点 $x_k=a+kh$，$h=\frac{b-a}{n}$，$k=0,1,...,n$。这里用 $T_n$ 表示用复化梯形法求得的积分值，其下表 n 表示等分数。

先考察一个字段 $[x_k,x_{k+1}]$，其中点 $x_{k+\frac{1}{2}}=\frac{1}{2}(x_k+x_{k+1})$，在该子段上二分前后的两个积分值

$$T_1=\frac{h}{2}[f(x_k)+f(x_{k+1})]$$
$$T_2=\frac{h}{4}[f(x_k)+2f(x_{k+\frac{1}{2}})+f(x_{k+1})]$$

显然有下列关系

$$T_2=\frac{1}{2}{T}_1+\frac{h}{2}f(x_{k+\frac{1}{2}})$$

将这一关系式关于 k 从 0 到 n - 1 累加求和，即可导出下列递推公式

$$T_{2n}=\frac{1}{2}{T}_n+\frac{h}{2}\sum _{k=0}^{n-1}f(x_{k+\frac{1}{2}})\text{\hspace{1em}(1)}$$

需要强调指出的是，上式中的 $h=\frac{b-a}{n}$ 代表二分前的步长，而

$$x_{k+\frac{1}{2}}=a+(k+\frac{1}{2})h$$

图 2-2 是变步长梯形法的算法框图。其中 $T_1$ 和 $T_2$ 分别代表二分前后的积分值。各框的含义是：

• [框 1] 准备初值
• [框 2] 按式（1）求二分后的梯形值

从第一个分点 $x=a+\frac{h}{2}$ 出发，取 h 为步长逐步向右跨，即可依次确定式（1）中的各个分点。图 2-2 将所得到的分点暂存于单元 x 中。

• [框 3] 控制精度
• [框 4] 修改步长

最终取二分后的积分值 $T_2$ 作为所求的结果。

例：用变步长方法计算

$$I={\int }_0^1\frac{\sin x}{x}dx$$

解：先对整个区间 [0, 1] 用梯形公式。$f(x)=\frac{\sin x}{x}$ 的数据参见下表。

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(1)=0.8414709$，故

$$T_1=\frac{1}{2}[f(0)+f(1)]=0.9207355$$

然后将区间二分，由于 $f(\frac{1}{2})=0.9588510$，利用递推公式（1）得

$$T_2=\frac{1}{2}{T}_1+\frac{1}{2}f(\frac{1}{2})=0.9397933$$

再二分一次，并计算新分点上的函数值

$$f(\frac{1}{4})=0.9896158$$

$$f(\frac{3}{4})=0.9088516$$

再用式（1）求得

$$T_4=\frac{1}{2}{T}_2+\frac{1}{4}[f(\frac{1}{4})+f(\frac{3}{4})]=0.9445135$$

这样不断二分下去，计算结果如下表所示（表中 k 代表二分次数，区间等分数 n = 2^k）。积分的准确值为 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;
}