Romberg算法

xiaoee posted @ 2009年4月01日 03:41 in Coding My Mind with tags 算法 数值计算 积分 递推 , 4524 阅读

Romberge算法是什么,其实是一种求数值积分的方法。就是先不断二分用递推梯形公式求出Tn,再通过Tn,T2n的线性组合求出复化Simpon积分Sn,再通过Sn,S2n的线性组合求出Cotes积分值,这样一步步往下算,它是一种加速算法。过程像下面这样:

                                                                                   T1 
                                                                                   T2    S1
                                                                                   T4    S2    C1
                                                                                   T8    S4    C2    R1
                                                                                   T16  S8    C4    R2   ....
 
下面是用C语言实现的代码:
#include"stdio.h"
#include"math.h"

#define A 0.//积分上界
#define B 1.//积分下界
#define EPS 0.00000001//精度

double results[12];

double F(double x){//被积函数
       return 4/(1+x*x);}

double T2n(double a,double b,int n){//梯形递推算法
       int i;
       double sum=0,h=(b-a)/pow(2,n),start=a+h;
       for(i=0;i<pow(2,n-1);i++)sum+=F(start+i*h*2);
       return sum*h+results[0]/2;
}

int Romberg(double a,float b,double eps){//Romberg算法
       int i,j,ok=0;
       double t,s;
       t=(b-a)*(F(a)+F(b))/2;
       results[0]=t;
       printf("%0.8f ",results[0]);
       for(i=1;!ok&&i<20;i++){
          t=T2n(a,b,i);
          printf("\n%0.8f ",t);
          for(j=0;!ok&&j<i;j++){
            s=(pow(2,2*(j+1))*t-results[j])/(pow(2,2*(j+1))-1);
            printf("%0.8f ",s);
            if(fabs(t-results[j])/(pow(2,2*(j+1))-1)<eps)ok=1;
            results[j]=t;
            t=s;}
            results[j]=t;}
       return j;         
}

int main(){
    printf("\n\nresult:%0.16f\n",results[Romberg(A,B,EPS)]);
    system("pause");
    return 0;
}
 

上面代码中是求4/(1+x^2从0到1的积分值,精确值是\pi

运行结果像下面的样子:


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter