mrom.zip.c

龙贝格算法的例子

# include<stdio.h>
# include<stdlib.h>
# include<math.h>
# define M 8
# define EPS 1.0e-6
double f(double x)
{return exp(x);}
double comp(double a,double b)
{static int first=0;
 static double Tn,h;
 double Hn,T2n,c=a;
 if(first==0)
 {h=b-a;
  Tn=0.5*h*(f(a)+f(b));
  first=1;
  return Tn;
  }
  Hn=0.0;h=h/2;
  while(c+h<b)
  {Hn=Hn+f(c+h); c=c+2*h;}
  Hn=Hn*h;
  T2n=0.5*Tn+Hn;
  Tn=T2n;
  return T2n;
}
int romberg(double a,double b,double T[M][M])
{int i,j,mi;
 double jj,err;
 for(i=0;i<M;i++)
 for(j=0;j<M;j++)
 T[i][j]=0.0;
 T[0][0]=comp(a,b);
 for(i=1;i<M;i++)
 {T[i][0]=comp(a,b);
   jj=1.0;
   for(j=1;j<=i;j++)
   {jj=jj*4.0;
     T[i][j]=(jj*T[i][j-1]-T[i-1][j-1])/(jj-1.0);
    }
    err=fabs(T[i][i-1]-T[i][i]);
    if(err<EPS){mi=i;goto endd;}
   }
   mi=M-1;
   endd:return(mi);
}
void main(void)
{int i,j,mi;
 static double a,b,T[M][M];
 printf("please input a:\n");
 scanf("%lf",&a);
 printf("please input b:\n");
 scanf("%lf",&b);
 mi=romberg(a,b,T);
 printf("\nRomberg method.\n");
  printf("\nintergral=%14.8f\n",T[mi][mi]);
  getch();
}