龙贝格求积.c

数值分析算法描述与习题解答

#include "math.h"
#include "stdio.h"
double romb(a,b,eps,f)
  double a,b,eps,(*f)();
{int m,n,i,k;
    double y[10],h,ep,p,x,s,q;
    h=b-a;
    y[0]=h*((*f)(a)+(*f)(b))/2.0;
    m=1; n=1; ep=eps+1.0;
    while ((ep>=eps)&&(m<=9))
      { p=0.0;
        for (i=0;i<=n-1;i++)
          { x=a+(i+0.5)*h;
            p=p+(*f)(x);
          }
        p=(y[0]+h*p)/2.0;
        s=1.0;
        for (k=1;k<=m;k++)
          { s=4.0*s;
            q=(s*p-y[k-1])/(s-1.0);
            y[k-1]=p; p=q;
          }
        ep=fabs(q-y[m-1]);
        m=m+1; y[m-1]=q; n=n+n; h=h/2.0;
      }
    return(q);
  }
  main()
  { double a,b,eps,t,rombf();
    a=0.0; b=0.35; eps=0.000001;
    t=romb(a,b,eps,rombf);
    printf("\n");
    printf("t=%e\n",t);
    printf("\n");
  }

  double rombf(x)
  double x;
  { double y;
    y=x/(x*x-4);
    return(y);
  }