几种数值积分方法

梯形多步法和辛普森积分
具体算法以及相关定理参考以下链接：

http://wenku.baidu.com/link?url=ZAScYogajXHrTTRa5xjpUPtS7OQQXZ_LfXaWNkczTtWf2MJgx0RZFUuca4iRGUcPHAVXOyGOIydwz3j-rriBiyLxRMnio7m4fNFcJOKJqui

#include <stdio.h>

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

#define EPS 1.0E-14        //计算精度
#define DIV 1000        //分割区间，值越大，精确值越高
#define ITERATE 20        //二分区间迭代次数，区间分割为2^n，初始化应该小一点，否则会溢出

#define RECTANGLE 0        //矩形近似
#define TRAPEZOID 1        //梯形近似
#define TRAPEZOID_FORMULA 2    //递推梯形公式
#define SIMPSON_FORMULA 3     //辛普森公式
#define BOOL_FORMULA 4        //布尔公式

double function1(double);       
double function2(double);
double function3(double);
void Integral(int, double f(double), double, double);      //矩形, 梯形逼近求定积分公式
void Trapezoid_Formula(double f(double), double a, double b);    //递推梯形公式
void Simpson_Formula(double f(double), double a, double b);        //递推辛普森公式
void Bool_Formula(double f(double), double a, double b);        //递推布尔公式
void Formula(int, double f(double), double, double);

double function1(double x)
{
    double y;
    y = cos(x);
    return y;
}

double function2(double x)
{
    double y;
    y=1/x;
   // y = 2+sin(2 * sqrt(x));
    return y;
}

double function3(double x)
{
    double y;
    y = exp(x);
    return y;
}

int main()
{    
    printf("y=cos(x) , x=[0, 1]\n");
    printf("Rectangle : "); Integral(RECTANGLE, function1, 0, 1);
    printf("Trapezoid : "); Integral(TRAPEZOID, function1, 0, 1);
    Formula(TRAPEZOID_FORMULA, function1, 0, 1);
    Formula(SIMPSON_FORMULA, function1, 0, 1);
    Formula(BOOL_FORMULA, function1, 0, 1);
    printf("\n");
    
    printf("y=1/x , x=[1, 5]\n");
    printf("Rectangle : "); Integral(RECTANGLE, function2, 1, 5);
    printf("Trapezoid : "); Integral(TRAPEZOID, function2, 1, 5);
    Formula(TRAPEZOID_FORMULA, function2, 1, 5);
    Formula(SIMPSON_FORMULA, function2, 1, 5);
    Formula(BOOL_FORMULA, function2, 1, 5);
    printf("\n");
    
    printf("y=exp(x) , x=[0, 1]\n");
    printf("Rectangle : "); Integral(RECTANGLE, function3, 0, 1);
    printf("Trapezoid : "); Integral(TRAPEZOID, function3, 0, 1);
    Formula(TRAPEZOID_FORMULA, function3, 0, 1);
    Formula(SIMPSON_FORMULA, function3, 0, 1);
    Formula(BOOL_FORMULA, function3, 0, 1);
    
    return 0;
}
//积分：分割-近似-求和-取极限
//利用黎曼和求解定积分
void Integral(int method, double f(double),double a,double b)   //f(double)，f(x), double a，float b，区间两点
{
    int i;
    double x;            
    double sum = 0;      //矩形总面积
    double h;            //矩形宽度
    double t = 0;        //单个矩形面积
    
    h = (b-a) / DIV;
     
    for(i=1; i <= DIV; i++)       
    {
        x = a + h * i;     
        switch(method)
        {
            case RECTANGLE :
                    t = f(x) * h;  
                break;
            case TRAPEZOID :
                    t = (f(x) + f(x - h)) * h / 2;
                break;
            default:
                break;
        }
        sum = sum + t;   //各个矩形面积相加
    }
    
    printf("the value of function is %.10f\n", sum);
}

//递推梯形公式
void Trapezoid_Formula(double f(double), double a, double b)//递推梯形公式
{
    unsigned int i, j, M, J=1;
    double h;
    double F_2k_1 = 0;
    double T[32];
    double res = 0;
    
    T[0] = (f(a) + f(b)) * (b-a)/2;
    for(i=0; i<ITERATE; i++)
    {
        F_2k_1 = 0;
        J = 1;
        
        for(j=0; j<=i; j++)    //区间划分
            J <<=  1;    //2^J
        h = (b - a) /J;    //步长
        //M = J/2;            //2M表示将积分区域划分的块数
        for(j=1; j <= J; j += 2)    //
        {
            F_2k_1 += f(a + j*h);    //f(2k-1)累加和
        }     
        T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
        res =  T[i+1];
        
        if(fabs(T[i+1] - T[i]) < EPS)
            if(i < 3) continue;
            else break;
    }
    
    printf("Trapezoid_Formula : the value of function is %.10f\n", res);
    //return res;
}
//辛普森公式
void Simpson_Formula(double f(double), double a, double b)        //辛普森公式
{
    unsigned int i, j, M, J=1;
    double h;
    double F_2k_1 = 0;
    double T[32], S[32];
    double res_T = 0, res_S = 0;
    
    T[0] = (f(a) + f(b)) * (b-a)/2;
    for(i=0; i<ITERATE; i++)
    {
        F_2k_1 = 0;
        J = 1;
        
        for(j=0; j<=i; j++)    //区间划分
            J <<=  1;    //2^J
        h = (b - a) /J;    //步长
        //M = J/2;            //2M表示将积分区域划分的块数
        for(j=1; j <= J; j += 2)    //
        {
            F_2k_1 += f(a + j*h);    //f(2k-1)累加和
        }     
        T[i+1] = T[i]/2 + h * F_2k_1;        //递推梯形公式
        res_T =  T[i+1];
        
        if(fabs(T[i+1] - T[i]) < EPS)
            if(i < 3) continue;
            else break;
    }
    
    for(i=1; i<ITERATE; i++)
    {
        S[i] = (4 * T[i] - T[i-1]) / 3;        //递推辛普森公式
        res_S = S[i];
        if(fabs(S[i] - S[i-1]) < EPS)
            if(i < 3) continue;
            else break;
    }
    
    printf("Simpson_Formula : the value of function is %.10f\n", res_S);
    //return res_S;
}

//布尔公式
void Bool_Formula(double f(double), double a, double b)        //布尔公式
{
    unsigned int i, j, M, J=1;
    double h;
    double F_2k_1 = 0;
    double T[32], S[32], B[32];
    double res_T = 0, res_S = 0, res_B;
    
    T[0] = (f(a) + f(b)) * (b-a)/2;
    for(i=0; i<ITERATE; i++)
    {
        F_2k_1 = 0;
        J = 1;
        
        for(j=0; j<=i; j++)    //区间划分
            J <<=  1;    //2^J
        h = (b - a) /J;    //步长
        //M = J/2;            //2M表示将积分区域划分的块数
        for(j=1; j <= J; j += 2)    //
        {
            F_2k_1 += f(a + j*h);    //f(2k-1)累加和
        }     
        T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
        res_T =  T[i+1];
        
        if(fabs(T[i+1] - T[i]) < EPS)
            if(i < 3) continue;
            else break;
    }
    
    for(i=1; i<ITERATE; i++)
    {
        S[i] = (4 * T[i] - T[i-1]) / 3;                //递推辛普森公式
        res_S = S[i];
        if(fabs(S[i] - S[i-1]) < EPS)
            if(i < 3) continue;
            else break;
    }
    
    for(i=1; i<ITERATE; i++)
    {
        B[i] = (16 * S[i] - S[i-1]) / 15;            //递推布尔公式
        res_B = B[i];
        if(fabs(B[i] - B[i-1]) < EPS)
            if(i < 3) continue;
            else break;
    }
    
    printf("Bool_Formula : the value of function is %.10f\n", res_B);
    //return res_B;
}

//采用二分区间的方法迭代，实际运行速度很慢
void Formula(int method, double f(double), double a, double b)//递推梯形公式
{
    unsigned int i, j, M, J=1;
    double h;
    double F_2k_1 = 0;
    double T[32], S[32], B[32];
    double res_B = 0, res_S = 0, res_T = 0;
    
    T[0] = (f(a) + f(b)) * (b-a)/2;
    for(i=0; i<ITERATE; i++)
    {
        F_2k_1 = 0;
        J = 1;
        
        for(j=0; j<=i; j++)    //区间划分
            J <<=  1;    //2^J
        h = (b - a) /J;    //步长
        //M = J/2;            //2M表示将积分区域划分的块数
        for(j=1; j <= J; j += 2)    //
        {
            F_2k_1 += f(a + j*h);    //f(2k-1)累加和
        }     
        T[i+1] = T[i]/2 + h * F_2k_1;        //递推公式
        res_T =  T[i+1];
        
        if(fabs(T[i+1] - T[i]) < EPS)
            if(i < 3) continue;
            else break;
    }
    
    switch(method)
    {
        default :
        case TRAPEZOID_FORMULA :
                printf("Trapezoid_Formula : the value of function is %.10f\n", res_T);
                //return res_T;
                break;
        case SIMPSON_FORMULA :
                for(i=1; i<ITERATE; i++)
                    {
                        S[i] = (4 * T[i] - T[i-1]) / 3;
                        res_S = S[i];
                        if(fabs(S[i] - S[i-1]) < EPS)
                            if(i < 3) continue;
                            else break;
                    }
                printf("Simpson_Formula : the value of function is %.10f\n", res_S);
                //return res_S;
                break;
        case BOOL_FORMULA :
                for(i=1; i<ITERATE; i++)
                {
                    S[i] = (4 * T[i] - T[i-1]) / 3;
                    res_S = S[i];
                    if(fabs(S[i] - S[i-1]) < EPS)
                        if(i < 3) continue;
                        else break;
                }
                for(i=1; i<ITERATE; i++)
                {
                    B[i] = (16 * S[i] - S[i-1]) / 15;
                    res_B = B[i];
                    if(fabs(B[i] - B[i-1]) < EPS)
                        if(i < 3) continue;
                        else break;
                }
                
                printf("Bool_Formula : the value of function is %.10f\n", res_B);
                //return res_B;
            break;
    }
    
}