数值积分实际上都是基于插值的，我们算的都是离散的一些个点与点对应的函数值，于是要求连续的积分是不可能的。

记等距分点a<x0<x1<...<b，距离h=(b-a)/n，n为分段数，设它的插值函数存在且为L(x)，我们用SL(x)dx的积分代替Sf(x)dx的积分值，这就是数值积分的最基本原理。

Newton-Cores公式，给出上面的积分为：

In=(b-a)∑[ C(n,k)f(xk) ]

其中C(n,k)是Cores系数，设x=a+th，则有

C(n,k)=(-1)^(n-k)/(nk!(n-k)! S∏(t-j)dt   j=0 to n,j!=k

当n=1时，算出C(1,0)=C(1,1)=1/2，即熟悉的梯形公式

I1=(b-a)(f(a)+f(b))/2

当n=2时，算出C(2,0)=1/6,C(2,1)=4/6,C(2,2)=1/6，这时就是Simpson公式，记得原来数学分析书上有提到过，它是按二次曲线插值的，有比较好的效率

I2=(b-a)(f(a)+4f((a+b)/2)+f(b))/6

这样一直算下去，会有更好的公式，且当n为偶数时，具有n+1次代数精度，但是当n>=8时，系数出现负项，将会不稳定，因此实际中只使用低次的Newton_Cores公式。

解决的办法有几种，一种是进行复化的求积公式，跟低次分段插值一样，将积分区间分成若干小段，每小段做低次积分，再求和；另一种办法就是用所谓的变步长积分。

其实所谓变步长法，就是我们熟知的递推法，在这里将看到递推法的妙处。我们拿到一个积分，和被积区间(a,b)，如何确定h呢？递推，先用较大的h试，如果不满意就将它细分，直到满意为止。

如果一开始只有一个分点，再将它们二分得四个，一直二分下去，那么原梯形公式将得到递推梯形公式，它是从n个分段二分到2n个的递推关系：

T2n=Tn/2+h/2∑f(xk+1/2)   [1]

右边求和部分是新插入的结点，于是就很方便程序控制，误差控制就考虑|T2n-Tn|的值，事后估计误差。但是直接用梯形法，收敛速度很慢，进一步优化发现，梯形法的余项展开为：

T(h)=I+a1h^2+a2h^4+a3h^6+...

如果h比较小，那T(h)=I+O(h^2)，为了提高效率，以h/2代h，得

T(h/2)=I+a1/4h^2+a2/16h^4+...

于是T1(h)=4/3T(h/2)-1/3T(h)=I+b1h^4+b2h^6+b3h^8+...，即T1(h)=I+O(h^4)，如果将前面公式代进来发现其实T1就是Simpson公式，同样的方法还可以继续。按这样的方法递推求积分通常称为Richardson外推加速法。设T(m,k)表示T(0,k)梯形值的m次加速值，那么有递推公式：

T(m,k)=4^m/(4^m-1) T(m-1,k+1) - 1/(4^m-1) T(m-1,k)   [2]

最后拿到手的公式就只有[1],[2]，按这样的方法构造一个T数表（它是一个三角形数表），得到的计算机上求积分的方法，就叫Romberg算法，算法简单描述为：
1，准备初值，T(0,0)=(b-a)(f(a)+f(b))/2
2，求梯形值，由公式[1]
3，求加速值，由公式[2]
4，精度控制

[T数表]
T(0,0)
T(0,1)T(1,0)
T(0,2)T(1,1)T(2,0)
T(0,3)T(1,2)T(2,1)T(3,0)
...

感觉这离最终实现还有一点点远的路，拿用C/C++实现吧，具体的实现方法还须考虑一些细节问题：
1，空间复杂度，T数表是一个三角形，如果开一个三角形大将是O(n*n)复杂度的，比较浪费，其实事实上有用的只是当前一行和下一行，于是只需要O(n)记录当前行，即一个梯形值和它的全部外推值。T[0]：当前梯形值，T[1]->T[m]外推值。
2、公式[2]两个系数比较复杂，如果直接算效率太低，而且中间结果可能膨胀，于是还需要算它的递推公式：
设T(m,k)=k1T(m-1,k+1)-k2T(m-1,k)
不难发现
k1/k2=4^m   (1)
k1=k2+1     (2)
(1)两边乘个4得到，4k1/k2=4^(m+1)，于是
k2'=1/(4k1/k2-1)=k2/(4k1-k2)
k1'=k2'+1
初值k1=4/3,k2=1/3

最后得到C/C++描述的实现程序：
#include <math.h>
#include <iostream.h>

double my_f(double x)
{
 return x==0.0?1.0:sin(x)/x;
}

double Romberg(double (*f)(double),double a,double b,double eps)
{
 double T[64];
 double h=b-a;
 int n=1;
 T[0]=h*(f(a)/4.0+f(b)/4.0+f(a+h/2.0)/2.0);//复化梯形公式
 T[1]=h*(f(a)/6.0+f(b)/6.0+f(a+h/2.0)/1.5);//辛甫生公式
 for(int i=2;fabs(T[i-2]-T[i-1])>eps;++i)
 {
  //计算递推梯形值,base
  h/=2.0;
  n<<=1;//分点数
  int j;
  double base=T[0]/h;
  double x=a+h/2.0;
  for(j=0;j<n;++j)
  {
   base+=f(x);
   x+=h;
  }
  base=base*h/2.0;
  //计算外推加速值,T[0]->T[i]
  double k1=4.0/3.0,k2=1.0/3.0;
  for(j=0;j<i;++j)
  {
   double hand=k1*base-k2*T[j];
   T[j]=base;
   base=hand;
   k2=k2/(4.0*k1-k2);
   k1=k2+1.0;
  }
  T[i]=base;
  //cout<<T[i]<<endl;
 }
 return T[i-1];
}

int main(void)
{
 cout<<Romberg(my_f,0.0,1.0,1e-8)<<endl;
 return 0;
}

如果积分区间比较大，当你把中间结果输出观察时，会发现有不稳定的现象，什么原因？由前面的梯形法的余项看，它是关于h的展开式，如果h的绝对值小于1，它才收敛，区间太大在一开始会有很大的h，自然结果是不稳定的，但是随着区间的二分，h终会达到收敛半径之内，从而使积分值收敛。

算法的复杂性，这个比较难描述，它和积分区间有关，如果限定只在[0,1]上积分，那设梯形值递推了n次，第i次将做2^i次的分点处值的计算，i次外推计算，最终计算次数为O(2^n)量级，但是分段半径h也会变成0.5^n，同时精度保证在h^(2n+4))的范围，即达到(0.5)^(2n^2+4n)，即达O(n*n)量级个二进制有效位，如果要达到十进制小数点后30位的有效数字的精度，因为100log<10,2>=30.x，故只需n=10左右即可，而做的运算次数只有2^10=1k次左右，可见效率很可观。

rickone 2006/10/19

评论