Sample solution in C for Homework 2, Qu, 2. Simpson's Rule.

======================================================================
Computer Code:
======================================================================
#include "math.h"
#include "stdio.h"

double f(double x)			// define the function to be integrated
{
    return exp(-x) * sin(x);
}

main()
{
    double x, a, b, h, simp, prevsimp, eps;
    int i, n, iter;

    eps = 0.000001;			   // set the precision; I've given 10^{-6}
    simp = 10000;
    b = 2.0;				   // set upper limit
    a = 0;				   // set lower limit

    printf ("         h          ans      ans - prev ans  (ans -prev ans)/h^4 \n");
    n = 2;
    while (fabs(prevsimp - simp) > eps)   // keep doubling no. of intervals until 
    {					  // desired precision has been obtained
	prevsimp = simp;
	h = (b - a) / n;		  // set width of one interval
	simp = f(a) + f(b) + 4*f(a+h);	  // initialize simp
	x = a + h;			  // set x equal to lower limit + h
	for (i = 2; i < n; i+=2)	  // sum over the remaining points
	{
	    x = x + h;			  // get current value of x
	    simp += 2 * f(x);		  // add an even point
	    x = x + h;			  // get current value of x
	    simp += 4 * f(x);		  // add an odd point
	}
	simp *= (h / 3);		  // multiply by h/3 (just once)
	if (n > 2)
	{printf ("%13.5f %13.8f %13.8f %13.8f \n", h, simp, (simp - prevsimp),
			(simp - prevsimp)/(h*h*h*h));} 
	else
	{printf ("%13.5f %13.8f \n", h, simp);}
	n *= 2;				  // double n for next iteration
    }
}
======================================================================
I keep doubling the number of intervals until the difference in the last
two results is less than the desired precision. I took 10^{-6} rather
than 10^{-4} asked for in the question; I also printed the results to 8
decimal places; note how this is done in the prinf statement.

I also printed out the difference in the last two results divided by
h^4. This is seen to go to a constant (as expected because Simpson is a
fourth order method). In fact the constant should equal

  15
- --- [f'''(2) - f'''(0)] 
  180

where 15 comes because we are taking the difference between the error for 2h
an the error for h, and the rest is the error expression quoted in
class.  This gives 0.155543... which seems to be consistent with the
numerics.

Overall, note the rapid convergence as h decreases.

======================================================================
Output:
======================================================================
         h          ans      ans - prev ans  (ans -prev ans)/h^4 
      1.00000    0.45376651 
      0.50000    0.46593497    0.01216846    0.19469530 
      0.25000    0.46658840    0.00065344    0.16727995 
      0.12500    0.46662712    0.00003872    0.15858173 
      0.06250    0.46662950    0.00000239    0.15630902 
      0.03125    0.46662965    0.00000015    0.15573505