变步长辛普森求数值积分

变步长辛普森求积算法，是计算数值积分里面的一个常用算法。与变步长梯形公式一样，将求积区间逐步二分，在每个子区间上运用辛普森公式，最后累加起来。与变步长梯形公式相比，该算法就是运算次数较少，比复合梯形公式计算速度快。但是精度是一样的。

下面贴代码：

/*变步长辛普森公式
*double simp(double a,double b,double eps,double (*f)(double))
*double a:区间的左端点
*double b:区间的右端点
*double eps:误差精度
*double (*f)(double)：积分函数
*返回值为double
*
*/

#include<iostream>
#include<cmath>
#include<iomanip>


double simp(double a, double b, double eps, double(*f)(double))
{
    double fa = (*f)(a);
    double fb = (*f)(b);
    double h = b - a;
    double T = h / 2.0 * (fa + fb);
    double x,S1 = T,S2;
    double T2;
    double ep = eps + 1.0;
    int n=1;
    while (ep >= eps)
    {
        double  p = 0.0;
        for (int k = 0; k <= n-1; k++)
        {
            x = a + (k + 0.5)*h;
            p = p + (*f)(x);
        }

        T2 = (T + h*p) / 2.0;
        S2 = (4.0*T2 - T) / 3.0;
        ep = fabs(S2 - S1);
        T = T2;
        S1 = S2;
        n = n * 2;
        h = h / 2.0;
      }
    return S2;
}


//测试函数: f(x)=exp(-x^2)

double func(double x)
{
    return exp(-1*pow(x,2));
}


//测试函数 f(x)=1/(1+x*x)
double func1(double x)
{
    return 1 / (1 + x*x);
}


int main()
{
    double a = 0.0;
    double b = 1.0;
    double eps = 0.00000001;
    double t = simp(a, b, eps, func1);
    std::cout << std::setprecision(8)<<"the value is " << t << std::endl;
    return 0;
}