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

谁践踏了优雅 2023-09-23 20:00 109阅读 0赞

![在这里插入图片描述][484e3ed8b1b540f481ccaee0b1d4ca43.png]

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;
        }
    }

运行结果得：  
![在这里插入图片描述][a29176e51ad94e2a9cb844769c82c7f7.png]

[484e3ed8b1b540f481ccaee0b1d4ca43.png]: https://img-blog.csdnimg.cn/484e3ed8b1b540f481ccaee0b1d4ca43.png
[a29176e51ad94e2a9cb844769c82c7f7.png]: https://img-blog.csdnimg.cn/a29176e51ad94e2a9cb844769c82c7f7.png