龙贝格python求积分 龙贝格算法求积分

---

梯形法 的算法简单，但精度低，收敛速度缓慢。如何提高收敛速度以节省计算量，自然是人们极为关心的课题。

根据梯形法的误差公式

$I-T_{n}\approx -\frac{h^{2}}{12}[f$

积分值 $T_{n}$ 的截断误差大致与 h² 成正比，因此步长减半后误差将减至 $\frac{1}{4}$，即有

$\frac{I-T_{2n}}{I-T_{n}}\approx \frac{1}{4}$

将上式移项整理，知

$I-T_{2n}\approx \frac{1}{3}(T_{2n}-T_{n})\tag1$

由此可见，只要二分前后两个积分值 $T_{n}$ 与 $T_{2n}$ 相当接近，就可以保证计算结果 $T_{2n}$

按式（1），积分值 $T_{2n}$ 的误差大致等于 $\frac{1}{3}(T_{2n}-T{n})$，如果用这个误差值作为 $T_{2n}$

$\overline T=T_{2n}+\frac{1}{3}(T_{2n}-T_{n})=\frac{4}{3}T_{2n}-\frac{1}{3}T_{n}\tag2$

应当是更好的结果。

再考察 此例，所求得的两个梯形值 $T_{4}=0.944 513 5$ 和 $T_{8}=0.945 690 9$

$\overline T=\frac{4}{3}T_{8}-\frac{1}{3}T_{4}=0.946 083 4$

却有六位有效数字。

按式（2）组合得到的近似值 $\overline T$，其实质究竟是什么呢？直接验证易知

$S_{n}=\frac{4}{3}T_{2n}-\frac{1}{3}T_{n}\tag3$

这就是说，用梯形法二分前后两个积分值 $T_{n}$ 与 $T_{2n}$ 按式（2）作线性组合，结果得到辛普生的积分值 $S_{n}$。

再考察辛普生法。按式

$I-S_{n}\approx -\frac{1}{180}(\frac{h}{2})^{4}[f$

其截断误差与 $h_{4}$

$\frac{I-S_{2n}}{I-S_{n}}\approx \frac{1}{16}$

由此得

$I\approx \frac{16}{15}S_{2n}-\frac{1}{15}S_{n}$

不难验证，上式右端的值其实等于 $C_{n}$，就是说，用辛普生法二分前后的两个积分值 $S_{n}$ 与 $S_{2n}$，按上式再作线性组合，结果得到柯特斯法的积分值 $C_{n}$，即有

$C_{n}\approx \frac{16}{15}S_{2n}-\frac{1}{15}S_{n}\tag4$

重复同样的手续，依据柯特斯的误差公式

$I-C_{n}\approx -\frac{2}{945}(\frac{h}{4})^{6}[f^{(5)}(b)-f^{(5)}(a)]$

可进一步导出下列龙贝格公式

$R_{n}=\frac{64}{63}C_{2n}-\frac{1}{63}C_{n}\tag5$

应当注意，龙贝格公式（5）已经不属于牛顿 - 柯特斯公式的范畴。

在步长二分的过程中运行公式（3）~（5）加工三次，就能将粗糙的积分值 $T_{n}$ 逐步加工成精度更高的龙贝格值 $R_{n}$，或者说，将收敛缓慢的梯形值序列 $T_{n}$ 加工成收敛迅速的龙贝格值序列 $R_{n}$，这种加速算法称为龙贝格算法，其加工流程如图 2-3 所示：

龙贝格算法的程序框图如图 2-4 所示：

例：基于 复化梯形的递推公式的变步长算法求积分 中的例子，用龙贝格算法计算

$I=\int_{0}^{1}\frac{\sin x}{x}dx$

得到的梯形值，计算结果如下表所示，表中 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;
}