Romberg 龙贝格算法求积分值

$Romberg$龙贝格算法

由误差的事后估计$$\frac{I-T_{2n}}{I-T_n}\approx \frac{1}{4}$$将上式移项得$$I-T_{2n}\approx \frac{1}{3}(T_{2n}-T_n)$$
也就是如果用$T_{2n}$作近似，则误差为约等于右边的部分，
那么将误差作为一种补偿加到$T_{2n}$上，那么得到的$$\overline{T}=\frac{4}{3}{T}_{2n}-\frac{1}{3}{T}_n$$应当是更好的结果。

经验证可知,实际上：
$$S_n=\frac{4}{3}{T}_{2n}-\frac{1}{3}{T}_n$$
即，补偿后的结果实际上是$Simpson$法的积分值，于是我们用梯形公式的递推可以直接算出辛普森的积分值。

相应的：$$Cotes积分值Cn\approx \frac{16}{15}{S}_{2n}-\frac{1}{15}{S}_n$$
同样的，依据$Cotes$法的误差公式，进一步导出$Romberg$公式
$$R_n=\frac{64}{63}{C}_{2n}-\frac{1}{63}{C}_n$$
由此，将收敛缓慢的梯形值序列加工成了收敛迅速的龙贝格值序列。

代码：

#include <iostream>
using namespace std;

double f(double x){
    if (x == 0) {
        return 1;
    }
    return sin(x)/x;
}

int main(int argc, const char * argv[]) {
    double b, a, epsilon;
    cout<<"请依次输入积分区间[a,b],允许的误差限epsilon"<<endl;
    cin>>a>>b>>epsilon;
    double h = b - a;
    double T1 = (f(a) + f(b)) * h / 2;
    double T2;
    double S1, S2;
    double C1, C2;
    double R1, R2;
    int k = 1; //k为二分次数
    double sum;
    double x;
    bool flg1 = false, flg2 = false, flg3 = false;
    do {
        if (flg3) {
            R1 = R2;
        }
        do {
            if (flg2) {
                C1 = C2;
            }
            do{
                if (flg1) {
                    ++k;
                    h /= 2;
                    T1 = T2;
                    S1 = S2;
                }
                sum = 0;
                x = a + h / 2;
                while (x < b) {
                    sum += f(x);
                    x += h;
                }
                T2 = sum * h / 2 + T1 / 2;
                S2 = T2 + (T2 - T1) / 3;
                flg1 = true;
            }while(k == 1);
            
            C2 = S2 + (S2 - S1) / 15;
            flg2 = true;
        } while (k == 2);
        R2 = C2 + (C2 - C1) / 63;
        flg3 = true;
    } while (k == 3 || fabs(R2 - R1) >= 3 * epsilon);
    
    cout<<R2<<endl;
}