romberg求积法终结版.c

来自「Romberg求积 Romberg求积」· C语言 代码 · 共 56 行

#include<stdio.h>
#include<math.h>
double f(double x)
{
 double y;
 y=4.0/(1.0+pow(x,2));
 return y;
}
main()
{
 double a;
 double b;
 double e;             //e为误差限
 int n;
 int i,k,m;
 double h;
 double F;
 double T[20][20];
 printf("请输入误差限ε:\n");
 scanf("%lf",&e);
 printf("请输入区间端点a,b(a<b):\n");
 scanf("%lf,%lf",&a,&b);
 h=(b-a)/2.0;
 T[0][0]=h*(f(a)+f(b));
 k=1;
 n=1;
 do
 {
  F=0;
  for(i=1;i<=n;i++)
  {
   F=F+f(a+(2*i-1)*h);
  }  //for
  T[0][k]=T[0][k-1]/2.0+h*F;
  for(m=1;m<=k;m++)
  T[m][k-m]=(pow(4,m)*T[m-1][k-m+1]-T[m-1][k-m])/(pow(4,m)-1);
  if(fabs(T[m-1][0]-T[m-2][0])<e)
  {
   printf("I≈%1.13lf.\n区间被划分为2*%d等份，每份的长度为%1.13lf.\n",T[m-1][0],n,h);
   break;
  }   //if
  else {
        h=h/2;
        n=2*n;
        k=k+1;
       }  //else
 }   //do
 while(1);
}   //main