算法原理

龙贝格（Romberg）求积法是甲酸定积分$T=\int_{a}^{b}f(x)dx$的最常用算法之一。
设$T_{m}(h)$ 步长为h时利用2m-2 阶牛顿-科茨（Newton-Cotes）公式计算得到的结果：$T_{m}(\frac{h}{2})$ 为将步长h减半后用2m-2 阶牛顿-科茨公式计算得到的结果。将他们机型线性组合，便得到步长为h的2m阶牛顿-科茨公式，即：

$$ T_{m+1}(h)=\frac{4^{m}T_{m}(h/2)-T_{m}(h)}{4^{m}-1} $$

其中$T_{1}(h)$为步长为h时的梯形公式计算得到的结果。$T_{1}(h/2)$为步长为h/2时的梯形公式计算得到的结果。
在实际计算时，龙贝格求积法按下表所示的计算格式进行，直到$|T_{m+1}(h)-T_{m}(h)| < \xi$为止。

算法实现

double CommonAlgorithm::integral(double (* vFunction)(double x), double a, double b, double eps = 1e-6){
    int m,n,i,k;
    double y[10],h,ep,p,x,s,q;

    // 迭代初值
    h=b-a;
    y[0]=h*(vFunction(a)+vFunction(b))/2.0;
    m=1; 
    n=1; 
    ep=eps+1.0;
    
    // 迭代计算
    while ((ep>=eps)&&(m<=9))
    { 
        p=0.0;
        
        for (i=0;i<=n-1;i++)
        { 
            x=a+(i+0.5)*h;
            p=p+vFunction(x);
        }

        p=(y[0]+h*p)/2.0;
        s=1.0;

        for (k=1;k<=m;k++)
        { 
            s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
        }

        ep=fabs(q-y[m-1]);
        m=m+1; 
        y[m-1]=q; 
        n=n+n; 
        h=h/2.0;
    }
    
    return(q);
}
 复制