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();
}