romberg.cpp

用龙贝格求积计算x*x*pow(M_E,-x*x)的积分值

#include <math.h>
#include <stdio.h>

#define eps 0.000001
#define max 20
#define M_E 2.71828182845904523536
double f(double x)
{
 if(x==0)
  return 1;
 else
  return (x*x*pow(M_E,-x*x));
}
void romberg(double a,double b)
{
 double t[max][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<max)&&(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++;
  h=h/2.0;
 }
 if(k>max)
  printf("method failed.\n");
 else
 {
  printf("\tt\t\ts\t\tc\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");
  }
 }
}

int main()
{
 double a=0.0,b=2.0;
 romberg(a,b);
 return 0;
}