龙贝格方法求解积分【Java实现】

public class Test {
    public static void main(String[] args) {
        double a, b, accuracy;   //a为区间下限，b为区间上限，accuracy为精度
        a =  0;//积分下界
        b =  0.8; //积分上界
        accuracy = 1e-8;//精度
        double h;   //步长
        h = b - a;
        double T1, T2;   //T1：二分前的梯形法积分值；T2：二分后的梯形法积分值；
        T1 = h / 2 * (1 + f(b));
        T2 = 0;
        int k = 1;  //记录二分次数
        System.out.println("T" +Math.pow(2, 0) +" = " +T1);
        double S1 = 0, S2;   //对T1与T2加权平均，得辛普森积分值
        double C1=0, C2;   //对S1与S2加权平均，得柯特斯积分值
        double R1 = 0, R2;   //对C1，C2加权平均，得龙贝格积分值
        int flag = 1;   //flag作为循环控制标志性变量
        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;    //计算梯形序列得下一个二分结果
            System.out.println("T" +Math.pow(2, k) + " = " + T2 );
            S2 = T2 + 1.0 / 3 * (T2 - T1);   //线性组合外推值simpson
            System.out.println("S"+Math.pow(2, k - 1) +" = " +S2);
            if (k == 1)   //至少外推2次得出S1，S2
            {
                k++;
                h /= 2;
                T1 = T2;
                S1 = S2;
                continue;
            }
            else
            {
                C2 = S2 + 1.0 / 15 * (S2 - S1);  //线性组合外推值Cotes
                System.out.println("C"+Math.pow(2, k - 2) +" = "+C2);

                if (k == 2)   //至少外推3次得出C1，C2
                {
                    C1 = C2;
                    k++;
                    h /= 2;
                    T1 = T2;
                    S1 = S2;
                    continue;
                }
                else
                {
                    R2 = C2 + 1.0 / 63 * (C2 - C1);   //线性组合外推至Romberg
                    System.out.println("R" +Math.pow(2, k - 3) + " = " + R2 );

                    if (k == 3)   //至少外推4次得出R1,R2
                    {
                        R1 = R2;
                        C1 = C2;
                        k++;
                        h /= 2;
                        T1 = T2;
                        S1 = S2;
                        continue;
                    }
                    else if (Math.abs(R2 - R1) >= accuracy)   //精度仍然不符合要求，继续二分步长、继续外推
                    {
                        R1 = R2;
                        C1 = C2;
                        k = k + 1;
                        h = h / 2;
                        T1 = T2;
                        S1 = S2;
                    }
                    else   //精度符合要求，修改flag为0，跳出while循环
                    {
                        flag = 0;
                        System.out.println("Romber算法求得数值积分结果为:" +R2 );
                    }
                }
            }
        }

    }
    static double f(double x)   //自定义被积函数
    {
        double result;
        result = Math.exp(Math.pow(-x,2));
        return result;
    }
}

运行结果得：