数字分析-数值积分-复化积分公式

复化积分公式应用计算以下公式，并能限制误差

$$\int_{0}^{1} \frac{sinx}{x} {\rm d}x$$

点击查看代码

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
/*************************************************
        积分公式数值实验
        1.复合梯形公式
        2.复合辛普森公式s
        0 - 1 区间 sinx / x 求积分

**************************************************/

/***************************************************
    函数功能：一重积分业务函数及其函数指针
    参数：x :未知数
    具体实现需自定义
****************************************************/
typedef double(*Fun)(double x);
double caculate_fun(double x)
{
    if(x == 0)
        return 1;
    else
        return (sin(x) / x);
}

/***************************************************
    函数功能：多重积分业务函数
    函数参数：   n:未知数个数
              x[n]:存储未知数数组
****************************************************/
double caculate_mfun(int n,double x[])
{
    int i;
    double f;
    f = 0.0;
    for(i=0; i<=n-1; i++)
    {
        f = f + x[i] * x[i];
    }
    return(f);
}
/**产生0-1之间均匀分布的随机数*/
double rnd1(double *r)
{
    int m;
    double s,u,v,p;
    s = 65536.0;
    u = 2053.0;
    v = 13849.0;
    m = (int)(*r/s);
    *r = *r-m*s;
    *r = u*(*r)+v;
    m = (int)(*r/s);
    *r = *r-m*s;
    p = *r/s;
    return(p);
}
/**************************************************
    自适应步长复合梯形公式
    误差控制计算公式|(T2n - Tn) / 3 <= jd|
    步长控制：不满足精度区间分段 n*2,h = 1/n
    Fun f为积分公式
    a-b 为积分区间
    n为分段数
    h为步长
    acc为累加项 accumulate累加和
    fa = f(a);
    fb = f(b);
    x_为累加项中的x
***************************************************/
double fuheTiXing(Fun f,double a,double b,double jd)
{
    double h,acc,fa,fb,x_,T;
    int i,n;
    double oldT;
    oldT = 0.0;  /*保存上次计算结果*/
    n = 2;       /*预设区间分两段如不满足要求继续分段*/
    do{
        h  = (b - a) / n; //步长
        fa = f(a);
        fb = f(b);
        acc = 0;
        x_  = a;
        T   = 0;
        for(i = 1;i < n;i++) //1 -> n-1  n-1次
        {
            x_  += h;
            acc += f(x_);
        }
        T = (h/2) * (fa + 2*acc + fb);
        if(fabs((T-oldT)/3) <= jd)  /*根据精度控制适应步长*/
        {
            break;              /*精度满足要求跳出循环*/
        }
        else{                   /*精度不满足要求分段*2继续计算一次*/
            oldT = T;
            n *=2;
        }
    }while(1);
    printf("复合梯形公式分段：%d ,满足精度要求\n",n);
    return T;
}
/************************************************************/
/**
    自适应步长复合辛普森求积公式
    匿名函数指针方便直接赋值函数名，但复杂函数函数参数会过长不利于查看
    |  |  |  |  |  |  |  |  |
    0  1  2  3  4  5  6  7  8
    x0    x1    x2    x3    x4
      h/2   h/2   h/2   h/2
    ——  —————  ————  —————        _x：从x0开始
    |     |     |     |     |
    0     2     4     6     8
    x0    x1    x2    x3    x4
          h
          —————  ————  —————      x_：从x1开始
    步长控制原理：根据误差要求及误差控制公式|(S2n - Sn) / 15| < jd
*/
/************************************************************/
double fuheSimpson(double (*f)(double),double a,double b,double jd)
{
    double h,fa,fb,acc,S,x_,_x;
    double oldS;
    int i,n;
    n = 2;       /*预设区间分两段如不满足要求继续分段*/
    oldS = 0.0;
    do{
        h = (b - a)/n;
        fa = f(a);
        fb = f(b);
        x_ = a;
        _x = a;
        acc = 0;
        for(i = 0;i < n;i++) //0 -> n-1共n次
        {
            if(i == 0)
            {
                _x  += 0.5*h;
                acc += 4*f(_x);
            }else
            {
                _x  += h;
                x_  += h;
                acc += 4*f(_x) + 2*f(x_);
            }
        }
        S = (h / 6)*(fa + acc + fb);
        if(fabs((S-oldS)/15) <= jd)  /*根据精度控制适应步长*/
        {
            break;                   /*精度满足要求跳出循环*/
        }
        else{                        /*精度不满足要求分段*2继续计算*/
            oldS = S;
            n *=2;
        }
    }while(1);
    printf("复合Simpson公式分段：%d ,满足精度要求\n",n);
    return S;
}
/****************************************************
    函数功能：蒙特卡洛多重积分计算函数
    参数说明：   n：积分计算的重数
              a[n]：存储积分下限值数组
              b[n]：存储积分上限值数组
                 f: 积分函数
*****************************************************/
double mtml(int n,double a[],double b[],double (*f)(int,double []))
{
    int m,i;
    double r,s,d,*x;
    x=malloc(n*sizeof(double));
    r=1.0;
    d=10000.0;
    s=0.0;
    for(m=0; m<=9999; m++)
    {
        for(i=0; i<=n-1; i++)
        {
            x[i] = a[i] + (b[i]-a[i]) * rnd1(&r);

        }
        s=s + (*f)(n,x)/d;
    }
    for(i=0; i<=n-1; i++)
    {
        s=s*(b[i]-a[i]);
    }
    free(x);
    return(s);
}
/*****************************************************
    函数功能：蒙特卡洛计算一重积分
    参数说明：a:积分下限
              b:积分上限
              f:积分函数
              rnd1：随机数函数(0-1)
*******************************************************/
double mtcl(double a,double b,double (*f)(double))
{
    int m;
    double r,d,x,s;
    r=1.0;
    s=0.0;
    d=10000.0;
    for(m=0; m<=9999; m++)
    {
        x = a + (b-a)*rnd1(&r);
        s = s + (*f)(x)/d;
    }
    s = s * (b-a);
    return(s);
}
int main(int argc, const char * argv[])
{
    double result_t,result_s,result_m;
    Fun fun = caculate_fun;
    result_t = result_s =0.0;
    result_m = 0.0;
    result_t = fuheTiXing(fun,0,1,0.000001);
    printf("复合梯形公式结果 = %f\n",result_t);
    result_s = fuheSimpson(caculate_fun,0,1,0.000001);
    printf("复合辛普森公式结果为：%f\n",result_s);
    result_m = mtcl(0,1,caculate_fun);
    //printf("蒙特卡洛积分结果为：%f\n",result_m);

    double a[3]={ 1.0,1.0,1.0};
    double b[3]={ 2.0,2.0,2.0};
    printf("\n");
   // printf("多重积分计算结果为s=%13.5e\n",mtml(3,a,b,caculate_mfun));
    return 0;
}