龙贝格求解数值积分

龙贝格求积公式，是一种梯形法的递推化，也是将求积区间逐次二分，将前一次二分的积分和Tn与下一次二分的积分和T2n有个递推关系，通过梯形公式的余项泰勒展开成级数形式，变量代换，进行整体的加减组合，能够得到关于T2n和Tn的线性组合S。同理，对S的线性组合可以推出C，对C的线性组合，可以推出R。精度在逐次提高，每一次递推提高2阶。

下面是实现代码：

/*龙贝格求积分
*double romb(double a,double b,double eps,double(*f)(double))
*double a:给出积分区间下限
*double b:积分区间上限
*double eps:误差精度
*double (*f)(double) :被积函数
*函数返回double
*/


#include<iostream>
#include<cmath>

double romb(double a, double b, double eps, double(*f)(double))
{
    int   i, k;
    double y[10], p, x,  q;
    double h = b - a;
    y[0] = h*((*f)(a) + (*f)(b)) / 2.0;
    int m = 1,n=1;
    double ep = eps + 1.0;

    while ((ep >= eps) && (m <= 9))
    {
        p = 0.0;
        for (i = 0; i < n; i++)
        {
            x = a + (i + 0.5)*h;
            p = p + (*f)(x);
        }
        p = (y[0] + h*p) / 2.0;
        for (k = 1; k <= m; k++)
        {
            q = (4.0*p - y[k - 1]) / 3.0;
            y[k - 1] = p;
            p = q;
        }
        ep = fabs(q - y[m - 1]);
        m = m + 1;
        y[m - 1] = q;
        n = n * 2;
        h = h / 2.0;
    }
    return q;
}


//积分函数:f(x)=x/(4+x*x)

double func(double x)
{
    return x / (4 + x*x);
}


int main()
{
    double a = 0.0;
    double b = 1.0;
    double eps = 0.0000001;
    double t = romb(a, b, eps, func);
    std::cout << "the value is " << t << std::endl;
    return 0;
}