算法原理
龙贝格（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$为止。
| 梯形法则 | 二阶公式 | 四阶公式 | 六阶公式 | 八阶公式 | ··· |
|---|---|---|---|---|---|
| $T_{1}(h)$ | |||||
| $T_{1}(\frac{h}{2})$ | $T_{2}(h)$ | ||||
| $T_{1}(\frac{h}{2^{2}})$ | $T_{2}(\frac{h}{2})$ | $T_{3}(h)$ | |||
| $T_{1}(\frac{h}{2^{3}})$ | $T_{2}(\frac{h}{2^{2}})$ | $T_{3}(\frac{h}{2})$ | $T_{4}(h)$ | ||
| $T_{1}(\frac{h}{2^{4}})$ | $T_{2}(\frac{h}{2^{3}})$ | $T_{3}(\frac{h}{2^{2}})$ | $T_{4}(\frac{h}{2})$ | $T_{5}(g)$ | ··· |
| : | : | : | : | : | : |
算法实现
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);
}