c龙贝格算法.txt

来自「C龙贝格算法 源程序 ｍａｔｌａｂ数值计算」· 文本 代码 · 共 58 行

龙贝格算 

--------------------------------------------------------------------------------

2007-03-06 22:06:00　标签：龙贝格算　　　[推送到技术圈] 


版权声明：原创作品，允许转载，转载时请务必以超链接形式标明文章 原始出处 、作者信息和本声明。否则将追究法律责任。http://lucky.blog.51cto.com/124662/19412 
#include<stdio.h>
#include<math.h>
  
double f(double x)
{
   return 4/(1+x*x);
}
 
double RombergInteg(f,a,b,eps)
   double (*f)();
   double a,b,eps;
{
   int i=1,j;
   double T[2][4],h=b-a;
   double x,s,r,det,m;
   T[0][0]=(f(a)+f(b))*h/2;
   while(1){
     s=0;
     for(x=a+h/2;x<b;x+=h) s+=f(x);
     T[1][0]=(T[0][0]+h*s)/2;
     m=4;
     for(j=1;j<=3&&j<=i;j++)
     {
       T[1][j]=(m*T[1][j-1]-T[0][j-1])/(m-1);
       m*=4;
     }
     if(i<=3)
     {
       r=T[1][i];
       det=fabs(r-T[1][i-1]);
     }
     else
     {
       r=T[1][3];
       det=fabs(r-T[0][3]);
     }
     if(det<eps) return r;
     for(j=0;j<4;j++) T[0][j]=T[1][j];
     h/=2; i++;
   }
}
 
main()
{
   double a=0,b=1,eps=1e-5,I;
   I=RombergInteg(f,a,b,eps);
   printf("I=%f\n",I);
}