【转】算定积分的程序

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double f (double x)
{
//    return sqrt(x); //f(x)=sqrt(x)

    return sin(x);
//    return x/(4+x*x);

}
//////////////////////////////////////

//梯形公式

double T(double x[],int n)
{
    return (x[n]-x[0])*(f(x[0])+f(x[n]))/2.0;
}

/////////////////////////////////////

//复化梯形公式

double ComT(double x[],double h,int n)
{
    int k;
    double m=0.0;
    for(k=1;k<n;k++)
        m+=f(x[k]);
    return h*(f(x[0])+2*m+f(x[n]))/2;
}

///////////////////////////////////

//辛普森公式

double S(double x[],int n)
{
    return (x[n]-x[0])*(f(x[0])+f(x[n])+4*f((x[0]+x[n])/2.0))/6.0;
}

////////////////////////////////////

//复化辛普森公式

double ComS(double x[],double h,int n)
{
    int k;
    double m=0.0,l=0.0;
    double *y;
    y= (double *)malloc(sizeof(double)*n);
    for(k=0;k<n;k++)
        y[k]=x[k]+h/2.0;
    for(k=0;k<n;k++)
        m+=f(y[k]);
    for(k=1;k<n;k++)
        l+=f(x[k]);
    return h*(f(x[0])+4*m+2*l+f(x[n]))/6.0;
    delete [] y;
}

///////////////////////////////////////////

//科特斯公式

double Cx(double x[],int n)
{
    double h;
    h=(x[n]-x[0])/4.0;
    double *y;
    y= (double *)malloc(sizeof(double)*5);
    for(int i=0;i<5;i++)
        y[i]=x[0]+i*h;
    return (y[4]-y[0])*(7*f(y[0])+32*f(y[1])+12*f(y[2])+32*f(y[3])+7*f(y[4]))/90.0;
//    delete [] y;

}

////////////////////////////////////////////

//复化科特斯公式

double ComCx (double x[],double h,int n)
{
    double m=0.0,t=0.0,l=0.0,s=0.0;
    int k;
    double *a,*b,*c;
    a=(double *)malloc(sizeof(double)*n);
    b=(double *)malloc(sizeof(double)*n);
    c=(double *)malloc(sizeof(double)*n);
    for(k=0;k<n;k++)
    {
        a[k]=x[k+1]-3*h/4;
        b[k]=x[k+1]-2*h/4;
        c[k]=x[k+1]-1*h/4;
    }
    for(k=0;k<n;k++)
    {
        m+=f(a[k]);
        t+=f(b[k]);
        l+=f(c[k]);
    }
    for(k=1;k<n;k++)
        s+=f(x[k]);
    return h/90.0*(7*f(x[0])+32*m+12*t+32*l+14*s+7*f(x[n]));
}

////////////////////////////////////////

//龙贝格公式

double Romberg(double x[],double h,int n)
{
    double *y;
    y=(double *)malloc(sizeof(double)*(2*n+1));
    for(int k=0;k<2*n+1;k++)
    {
        if(k%2==0)
            y[k]=x[k/2];
        else
            y[k]=x[(k-1)/2]+h/2.0;
    }
    return 64/63*ComCx(y,h/2.0,2*n)-1/63*ComCx(x,h,n);
}

//////////////////////////////////////////////

//高斯求积公式

double Gauss(double x[],int n)
{
    return (x[n]-x[0])/2.0*(5/9.0*f((-1)*sqrt(15)*(x[n]-x[0])/10.0+(x[n]+x[0])/2.0)
         +8*f((x[n]+x[0])/2.0)/9.0+
         5/9.0*f(sqrt(15)*(x[n]-x[0])/10.0+(x[n]+x[0])/2.0));
}

int main()
{
    printf("请输入积分区间[a,b]:");
    double a,b;
    scanf("%lf%lf",&a,&b);
    printf("请输入区间数n:");
    int n;
    scanf("%d",&n);
    double h;
    h=(b-a)/(1.0*n);
    printf("区间的步长为:%lf\n",h);
    double *x;
    x=(double *)malloc(sizeof(double)*(n+1));
    for(int i=0;i<=n;i++)
        x[i]=a+i*h;
    double r;
    r=T(x,n);
    printf("梯形公式的结果是:%lf\n",r);
    printf("\n");
    r=ComT(x,h,n);
    printf("复化梯形公式结果是:%lf\n",r);
    printf("\n");
    r=S(x,n);
    printf("辛普森公式的结果是:%lf\n",r);
    printf("\n");
    r=ComS(x,h,n);
    printf("复化辛普森公式的结果是:%lf\n",r);
    printf("\n");
    r=Cx(x,n);
    printf("科特斯公式的结果是:%lf\n",r);
    printf("\n");
    r=ComCx(x,h,n);
    printf("复化科特斯公式的结果是:%lf\n",r);
    printf("\n");
    r=Romberg(x,h,n);
    printf("龙贝格公式的结果是:%lf\n",r);
    printf("\n");
    r=Gauss(x,n);
    printf("高斯公式的结果是:%lf\n",r);
    printf("\n");
    
    system("pause");
    return 0;
}