romberg2.c

(1) 设计算法并编制程序

#include <math.h>
#include <stdio.h>
#define eps 0.00000001
#define N  30
#define   PI 3.1415926
double f(double x)
{ double y=0.9;
return  sqrt(1-y*y*cos(x)*cos(x)) ;
}
void romberg(double a,double b)
{
 double t[N][N]={0},h=b-a,e=1.0+eps;
 double new;
 int i,j,k=1,m;
 t[1][1]=h*(f(a)+f(b))/2.0;
 while((k<N)&&(e>eps))
{new=0;
  for(i=1;i<=(int)(pow(2,k-1));i++)
   new=new+f(a+(i-0.5)*h);
  t[k+1][1]=(t[k][1]+h*new)/2.0;
  for(m=1;m<=k;m++)
  {
   if(m>3)
    break;
   t[k+1][m+1]=(pow(4,m)*t[k+1][m]-t[k][m])/(pow(4,m)-1);
  }
  if(k>=4)
   e=fabs(t[k+1][3]-t[k][3]);
  k++;
  h=h/2.0;
 }
 if(k>N)
  printf("method failed.\n");
 else
 {
  printf("\tt(0)(k)\t\tt(1)(k-1)\tt(2)(k-2)\tt(3(k-3)\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+1][j+1]);
    printf("\n");
  }
 }
printf("T=%0.9lf\n",t[k][m]);
}
int  main()
{
 double a=0,b=2*PI;
romberg(a,b);
 return 0;  }