龙贝格求积分算法C语言实现

数值分析中的龙贝格算法，是用来求积分的，算法主要是通过不停的迭代，使得最终积分函数值和真实函数值逼近。

对于足够“光滑”的函数来说，收敛速度比较快

 

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define M_PI       3.14159265358979323846
#define M_E        2.71828182845904523536
#define FUN(x)        pow(4 - pow(sin((x)),2),0.5)
#define FUN2(x)     ((x)*pow((x)*(x) + 1,0.5))//积分函数
#define N 20  //最大加速次数
#define A 0  //积分区间
#define B 3//积分区间
#define FAC        1//积分函数前面的系数
#define T_ZERO  1.423024947075771e+001//第一次用梯形公式时积分函数值
void main()
{ 
    double T[N][N];
    int k,m;
    int i,j;
    T[0][0] = T_ZERO;
    T[1][0] = 6.866658294205204e-001;
    for(k = 1;k<N;k++)
    {
        double t = 0,h = (B-A) * pow(0.5,k);
        int n = (int)pow(2,k-1);
#ifdef DEBUG_M
        printf("%d\t %lf\n",n,h);
#endif
        for( m = 0;m < n;m++)
        {
            t += h  * FUN2(h + 2*h*m);
        }
        T[k][0] = 0.5*T[k-1][0] + t;
    }
    for(k = 1;k<N;k++)
    {
        for(m = 1;m<=k;m++)
            T[k][m] = pow(4,m)/(pow(4,m)-1) *  T[k][m-1] -  1/(pow(4,m)-1) * T[k-1][m-1];
        if(fabs(T[k][k] - T[k-1][k-1]) < 1e-5)
            break;
    }
    for(i = 0;i<=k;i++)
        for(j = 0;j<=i;j++)
            T[i][j] *= FAC;
    for(i = 0;i<=k;i++)
    {
        for(j = 0;j<=i;j++)
            printf("%.15lf\t\t",T[i][j]);
        printf("\n");
    }
        
}