PHYSICS 115/242

    Example solution in C, for Qu.1, HW 2

===========================================================================
Source Code
===========================================================================
#include "math.h"
#include "stdio.h"

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

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

    TWO = 2;				   // initialization
    exact = log(TWO);			   // exact answer
    ITERMAX = 8;			   // set no. of iterations
    b = 1.0;				   // set upper limit
    a = 0;				   // set lower limit

    printf ("         h          ans      (exact-ans)/h^2 \n");
    n = 1;
    for (iter = 0; iter < ITERMAX; iter++) // keep doubling no. of intervals
    {
	h = (b - a) / n;		  // set width of one interval
	trap = 0.5 * (f(a) + f(b));	  // initialize trap
	x = a;				  // set x equal to lower limit
	for (i = 1; i < n; i++)		  // sum over n intervals
	{
	    x = x + h;			  // get current value of x
	    trap += f(x);		  // add f(x) to trap
	}
	trap *= h;			  // multiply by h (just once)
	n *= 2;				  // double n for next iteration
	printf ("%12.5f %12.5f %12.5f \n", h, trap, (trap - exact)/(h*h));
    }
}

===========================================================================
Output:
===========================================================================

         h          ans      (exact-ans)/h^2 
     1.00000      0.75000      0.05685 
     0.50000      0.70833      0.06074 
     0.25000      0.69702      0.06203 
     0.12500      0.69412      0.06238 
     0.06250      0.69339      0.06247 
     0.03125      0.69321      0.06249 
     0.01562      0.69316      0.06250 
     0.00781      0.69315      0.06250 

===========================================================================
Comments:
===========================================================================

One sees from the last column that, for small h, the error is 0.0625 h^2
= h^2 / 16. This agrees with the error estimate given in the question
since f'(x) = -1/(1+x)^2.