romberg.c

计算方法中romberg算法的c程序源码

#include <math.h>
double f();
double Romberg();
main()
{
 double a,b,EPS;
 int FAZ;
 printf("\nInput a,b,EPS,FAZ\n");
 scanf("%lf,%lf,%lf,%d",&a,&b,&EPS,&FAZ);
 Romberg(a,b,EPS,FAZ);
}

double f(double x)
{
 return(sin(x));
}

double Romberg(a,b,EPS,FAZ)
double a,b,EPS;
int FAZ;
{
 int FAG,i=1,j,k,m;
 double y1,y2,H,New,F,ER,x,s;
 double S[2][16];
 printf("FAZ  s              ER               FAG\n");
 y1=f(a); 
 y2=f(b);
 S[0][1]=(b-a)*(y1+y2)/2;
 H=b-a;
 while(1)
  {
   New=0;
   for(j=1;j<=pow(2,i-1);j++)
    {
     x=a+(2*j-1)*H/2;
     F=f(x);
     New=New+F;
    }
   S[1][1]=(S[0][1]+H*New)/2;
   for(k=1;k<=i;k++)
    S[1][k+1]=(pow(4,k)*S[1][k]-S[0][k])/(pow(4,k)-1);
   ER=fabs(S[1][i+1]-S[0][i]);
   if(ER>EPS && i<=FAZ)
    {
     for(m=1;m<=i+1;m++)
      S[0][m]=S[1][m];
     i=i+1;
     H=H/2;
    }
   else
    {
     s=S[1][i+1];
     break;
    }
   }
   if(ER<=EPS)  FAG=1;
   else         FAG=-1;  
   printf("%d    %11.10lf   %11.10lf     %d",FAZ,s,ER,FAG);
   printf("\nReturn:press Space     Quit:press else\n");
   if(getch()==32)  return main();
}