[转贴]Romberg算法求数值积分的原理与实现

风花雪月

数值积分实际上都是基于插值的，我们算的都是离散的一些个点与点对应的函数值，于是要求连续的积分是不可能的。
记等距分点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++描述的实现程序：

1. #include <math.h>
   
2. #include <iostream.h>
   
3. double my_f(double x)
   
4. {
   
5. return x==0.0?1.0:sin(x)/x;
   
6. }
   
7. double Romberg(double (*f)(double),double a,double b,double eps)
   
8. {
   
9. double T[64];
   
10. double h=b-a;
    
11. int n=1;
    
12. T[0]=h*(f(a)/4.0+f(b)/4.0+f(a+h/2.0)/2.0);//复化梯形公式
    
13. T[1]=h*(f(a)/6.0+f(b)/6.0+f(a+h/2.0)/1.5);//辛甫生公式
    
14. for(int i=2;fabs(T[i-2]-T[i-1])>eps;++i)
    
15. {
    
16.   //计算递推梯形值,base
    
17.   h/=2.0;
    
18.   n<<=1;//分点数
    
19.   int j;
    
20.   double base=T[0]/h;
    
21.   double x=a+h/2.0;
    
22.   for(j=0;j<n;++j)
    
23.   {
    
24.    base+=f(x);
    
25.    x+=h;
    
26.   }
    
27.   base=base*h/2.0;
    
28.   //计算外推加速值,T[0]->T[i]
    
29.   double k1=4.0/3.0,k2=1.0/3.0;
    
30.   for(j=0;j<i;++j)
    
31.   {
    
32.    double hand=k1*base-k2*T[j];
    
33.    T[j]=base;
    
34.    base=hand;
    
35.    k2=k2/(4.0*k1-k2);
    
36.    k1=k2+1.0;
    
37.   }
    
38.   T[i]=base;
    
39.   //cout<<T[i]<<endl;
    
40. }
    
41. return T[i-1];
    
42. }
    
43. int main(void)
    
44. {
    
45. cout<<Romberg(my_f,0.0,1.0,1e-8)<<endl;
    
46. return 0;
    
47. }

复制代码

如果积分区间比较大，当你把中间结果输出观察时，会发现有不稳定的现象，什么原因？由前面的梯形法的余项看，它是关于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


转自算法驿站：http://blog.programfan.com/