romberg.c

数值分析中RomBerg插值算法

//用ROMBERG方法计算sinx/x在0-1上的积分
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#define MAX 20
void main() 
 { 
   double a=0,b=1,r=0.0000005;
   double h=(b-a)/2; 
   double T[MAX][MAX]; 
   int k=1,n=1,i,j,m; 
   double F=0.0; 
   double result=0.0;
   
   for(i=0;i<MAX;i++)
   {
     for (j=0;j<MAX;j++)
     {
       T[i][j]=0.0;
     }
   }
   T[0][0]=h*(1+sin(b)/b);
 do{
	 F=0;
   for (i=1;i<=n;i++) 
   { 
     F=F+sin(a+(2*i-1)*h)/(a+(2*i-1)*h); 
   } 
   
     T[0][k]=T[0][k-1]/2+h*F;
   
   for(m=1;m<=k;m++)
   {
     T[m][k-m]=(pow(4,m)*T[m-1][k-m+1]-T[m-1][k-m])/(pow(4,m)-1);
   }
   m=m-1;
   if(fabs(T[m][0]-T[m-1][0])>=r)
   {
		h=h/2;
		n=2*n;
		k=k+1;
   }
   else break;
   }
   while(fabs(T[m][0]-T[m-1][0])>=r);
   result=T[m][0];
   printf("%f",result); 
}