/******************************************************************

Computes numerically, using the midpoint rule, the integral of
1/(1 + x)^2 from 0 to 1. The exact result is 1/2.

The results, given at the bottom, converge rapidly to this as the
number of intervals, n, is increased. I keep doubling the number
intervals, so n = 1, 2, 4, 8, ...,  128.

The convergence to the exact answer is clearly seen,

******************************************************************/

#include <math.h>
#include <stdio.h>

double f(double x)			// define the function to be integrated
{
    return 1 / ((1+x) * (1+x));
}

main()
{
    double x, a, b, h, midpt, exact, error;
    int i, n, iter, ITERMAX;

    ITERMAX = 8;			   // set no. of iterations
    b = 1.0;				   // set upper limit
    a = 0;				   // set lower limit
    exact = 0.5;                           // exact answer

    printf ("      n        h            ans         error \n");
    n = 1;
    for (iter = 0; iter < ITERMAX; iter++) // keep doubling no. of intervals
    {
	h = (b - a) / n;		  // set width of one interval
	midpt = 0 ;	  		  // initialize midpt
	x = a + 0.5 * h;		  // set x equal to lower limit + (1/2)h
	for (i = 1; i <= n; i++)	  // sum over n intervals
	{
	    midpt += f(x);		  // add f(x) to midpt
	    x += h;			  // add h to x for next evaluation
	}
	midpt *= h;			  // multiply by h (just once)
	error = midpt - exact;
	printf (" %6d  %12.6f %12.6f %12.6f \n", n, h, midpt, error);
	n *= 2;				  // double n for next iteration
    }
}
/*
            OUTPUT
      n        h            ans         error 
      1      1.000000     0.444444    -0.055556 
      2      0.500000     0.483265    -0.016735 
      4      0.250000     0.495548    -0.004452 
      8      0.125000     0.498867    -0.001133 
     16      0.062500     0.499716    -0.000284 
     32      0.031250     0.499929    -0.000071 
     64      0.015625     0.499982    -0.000018 
    128      0.007812     0.499996    -0.000004 
*/