Listing 1

#include <math.h>

/* Adaptive quadrature routine.
   Called with:
        1ft:          left endpoint
        rt:           right endpoint
        f:            pointer to function to integrate
        err:          pointer to allowable error

   Returns
        return value  evaluated integral
        err           pointer to upper limit of actual error
*/

double adaptive (double lft, double rt, double (*f)(double x), double *err)
{
   double h;                         /* Width of this interval. */
   double p1, p2, p3, p4, p5;        /* Value of function at the 5 points. */
   double error;                     /* Actual error on this interval. */
   double err1, err2;                /* Errors of the sub-intervals. */
   double int3,                      /* Value of the integral with Simpson */
         int5;                      /* 3- and 5-point approximations. */

   h = rt - lft;                     /* Get the width of the interval. */
                        /* Get the 2 approximations to the integral.*/

   p1 = (*f) (lft);
   p2 = (*f)(lft + h/4.0);
   p3 = (*f)(lft + h/2.0);
   p4 = (*f)(rt - h/4.0);
   p5 = (*f)(rt);

   int3 = h * (p1 + 4.0 * p3 + p5) / 6.0;
   int5 = h * (p1 + 4.0 * p2 + 2.0 * p3 + 4.0 * p4 + p5) / 12.0;

   error = fabs (int3 - int5) / 15.0;  /* Find the actual error. */
   if (error < *err)                   /* If within tolerance, */
   {
      *err = error;                   /* copy the actual error to caller, */
      return (int5);                  /* return the more accurate approx. */
   }
   else                                /* Otherwise, */
   {
      err1 = err2 = *err / 2.0;       /* cut the allowable error in half, */
      int5 = adaptive (lft, lft + h/2.0, f, &err1);  /* integrate the two */
      int5 += adaptive (lft + h/2.0, rt, f, &err2);  /* half-regions, */
      *err - err1 + err2;             /* sum the real errors, */
      return (int5);                  /* and return the integral. */
   }
}

/* End of File */