数值计算（三）之数值积分（变步长复化梯形公式，龙贝格公式）

数值积分

有一些复杂的函数，或者一些简单的函数但找不到原函数，这种情况下，我们就可以使用一些方法近似求得积分值。

这里三种近似的求积公式

这里介绍变步长复化梯形算法和龙贝格算法

变步长复化梯形算法

基本原理就是在求积区间内用梯形公式求积分，如果精度不够，我们就把区间长度对半分，更加细致，重复计算。直至本次计算的结果与之前计算的结果之差满足精度。

变步长复化梯形公式

double VariableStepSize(double (*f)(double),double a,double b,double e)
{//f为求积函数,a和b为求积区间,e是精确要求
    long long i,times=1;//用来记录增加节点的个数
    double ans,t,sum;
    ans=(f(a)+f(b))*(b-a)/2;//利用梯形公式粗略计算积分
    double numerator,denominator=1;
    do//至少循环一次
    {
        t=ans;
//        printf("t=%.7lf\n",t);
        sum=0;
        numerator=1;    denominator*=2;
        for(i=0;i<times;i++)
        {//新增加的点用了两次
            sum=sum+2*f( numerator/denominator );
            numerator+=2;
        }//End for-i
        times*=2;
        sum=(sum)/((b-a)*times*2);//新增加的部分
        ans=t/2+sum;
    }while((ans-t-e)>=0);//满足精度要求以后退出
    return ans;
}

龙贝格算法

龙贝格算法则是利用了梯形公式，辛普森公式，柯特斯公式之间的关系，逐步增加精度

龙贝格算法

double Longberg(double (*f)(double),double a,double b,double e)
{//f为求积函数,a和b为求积区间,e是精确要求
    int i,j,k=1;//每次递推都可以推导一行
    int times=1;
    double sum;
    double table[10][10]={0};//递推的二维矩阵,简单地固定位10*10
    double coef[10][2]={0};//用来保存系数,比如说4/3,16/15
    double LastLeft,LastRight;//记录推导一行中最后两个数，用以判断精度
    double numerator,denominator=1;//分子分母
    table[0][0]=( f(a)+f(b)*(b-a) )/2;    numerator=4;
    for(i=1;i<10;i++)
    {
        coef[i][0]=numerator/(numerator-1);//前一个系数，比如说4/3
        coef[i][1]=(-1)/(numerator-1);//后一个系数,比如说-1/3
        numerator*=4;
    }//End for-i
    do
    {
        sum=0;
        numerator=1;    denominator*=2;
        for(i=0;i<times;i++)
        {
            sum=sum+f( numerator/denominator );
            numerator+=2;
        }
        times*=2;   sum=(sum)/((b-a)*times*2);
        table[k][0]=table[k-1][0]/2+sum;//求得Tn

        for(i=1;i<=k;i++)//向外递推过程
            table[k][i]=coef[i][0]*table[k][i-1]+coef[i][1]*table[k-1][i-1];
        LastLeft=table[k][k-1]; LastRight=table[k][k];
        k++;
    }while((LastRight-LastLeft-e)>=0);//满足精度则退出
    return LastRight;
}