龙贝格求积公式 c语言实现 数值积分

龙贝格求积公式 c语言实现 数值积分

标签：计算方法实验

---

/*
    本实验用龙贝格求积公式求4 / (1 + x^2)在[0, 1]的定积分。
*/
#include <stdio.h>
#include <math.h>

double f(double x){
    return 4 / (1 + x * x);
}

int main(){
    double a = 0, b = 1, eps = 0.00001;  //上下限，精度
    double s, x, s1, s2, c1, c2, r1, r2, h = b - a, t1 = h * (f(a) + f(b)) / 2.0, t2, temp;
    int k = 1;

    do{
        s = 0, x = a + h / 2.0;
        do{
            s += f(x);
            x += h;
        }while(x < b);
        t2 = (t1 + h * s) / 2.0;
        s2 = t2 + (t2 - t1) / 3.0;
        if(k == 1)  {t1 = t2, s1 = s2, h /= 2.0, k += 1;  continue;}
        c2 = s2 + (s2 - s1) / 15.0;
        if(k == 2)  {t1 = t2, s1 = s2, c1 = c2, h /= 2.0, k += 1;  continue;}
        r2 = c2 + (c2 - c1) / 63.0;
        if(k == 3)  {t1 = t2, s1 = s2, c1 = c2, r1 = r2, h /= 2.0, k += 1;  continue;}
        temp = r1, r1 = r2, t1 = t2, s1 = s2, c1 = c2, h /= 2.0, k += 1;
    }while(fabs(r1 - temp) >= eps);  //即fabs(r2 - r1) >= eps
    printf("answer = %f\n", r2);
    return 0;
}

实验结果