longbeigejifen.cpp

用c++编写的龙贝格积分的程序 附有实例和计算结果

#include<maTh.h>
#include<iosTream.h>
#include<sTdio.h>
#define c 20
#define EPS 0.0000001
      double  f(double x)
	  {
		  if(x==0) return (1);
		  else 
			  return (sin(x)/x);
	  }
	  void  romberg(double a, double b)
	  {
		  double T[c][4]={0},h=1.0, e=1.0+EPS;
		  double fnew;
		  int i,j,k=1,m;
		  T[0][0]=h*(f(a)+f(b))/2.0;
		  while((k<c)&&(e>EPS))
		  {	  fnew=0;
		  for(i=1; i<=(int)(pow(2,k-1));i++)
			  fnew=fnew+f(a+(i-0.5)*h);
		  T[k][0]=(T[k-1][0]+h*fnew)/2.0;
		  for(m=1;m<=k;m++)
		  {
			  if(m>3)break;
			  T[k][m]=(pow(4,m)*T[k][m-1]-T[k-1][m-1])/(pow(4,m)-1);

		  }
		  if(k>=4)
			  e=fabs(T[k][3]-T[k-1][3]);
		  k=k+1;
		  h=h/2.0;
	  }
	  if(k>c)
		  printf("Method Failed!\n");
	  else
	  {

		  printf("T\t\tS\t\t\C\t\tR\n");
		  for(i=0;i<k;i++)
		  {
			  printf("k=%d\t",i);
			  for(j=0;j<4;j++)
				  if(i>=j)
					  printf("%0.9lf\t",T[i][j]);
				  printf("\n");
		  }
	  }
	  }
			 
      void main()
	  {
		  cout<<"利用龙贝格积分法计算f(x)=sinx/x dx"<<endl;
			

           double a=0.0, b=1.0;
		   cout<<"在区间"<<"["<<a<<","<<b<<"]"<<"求出龙贝格积分"<<endl;
			   romberg(a,  b);
	  }