Romberg

Romberg Integration (statement only)

#include<stdio.h>
#include<conio.h>
#include<math.h>
#define F(x) 1.0/x
int main (void)
{
	int i,j,k,m,n;
   float a,b,h,sum,x,r[10][10];
   printf("End points of the interval:\n");
   scanf("%f%f",&a,&b);
   printf("Maximum number of times the sub interval bisected: ");
   scanf("%d",&n);
   /* compute area :: trapezoidal */
   h=b-a;
   r[1][1]=h*(F(a)+F(b))/2.0;
   printf("\n%f\n",r[1][1]);
   /* Start ROMBERG integration */
   for(i=2;i<=n+1;i++)
   {
   	m=pow(2,(i-2));
      h=h/2;
   /* use RECURSIVE trpezoidal rule */
      sum=0.0;
      for(k=1;k<=m;k++)
      {
      	x=a+(2*k-1)*h;
         sum=sum+F(x);
      }
      r[i][1]=r[i-1][1]/2.0+h*sum;
   /* Compute Richardson improvement */
      for(j=2;j<=i;j++)
      {
      	r[i][j]=r[i][j-1]+(r[i][j-1]-r[i-1][j-1])/(pow(4,j-1)-1);
      }
   /* Write result of improvment for ith refinement */
      for(j=1;j<=i;j++)
      	printf("\n%f",r[i][j]);
      printf("\n");
   /* Test for accuracy */
      if(fabs(r[i-1][i-1]-r[i][i])<0.000001)
   /* Stop further refinement */
      {
      	printf("\n ROMBERG INTEGRATION: %f\n",r[i][i]);
         goto stop;
      }
      else
      {
      	continue;
      }
   }
   /* Write final RESULT */
   printf("\n\tROMBERG INTEGRATION = %f \n",r[n+1][n+1]);
   stop:
   printf("\tEnd!");
   getch();
}