Listing 2

#include <math.h>

/* Adaptive quadrature routine.
   Called with:
        lft:           left endpoint
        rt:            right endpoint
        f:             pointer to function to integrate
        err:           pointer to allowable error
        p1:            function value at left endpoint
        p2:            function value at center of interval
        p3:            function value at right endpoint

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

static double adaptive2 (double lft, double rt, double (*f)(double x), double *err,
                         double p1, double p3, double p5)
{
   double p2, p4;                    /* Function values at pts 2 & 4. */
   double h;                         /* Width of this interval. */
   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. */
   p2 = (*f)(lft + h/4.0);
   p4 = (*f)(rt - h/4.0);

   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;       /* integrate the two half-regions, */
      int5 = adaptive2 (lft, lft + h/2.0, f, &err1, p1, p2, p3);
      int5 += adaptive2 (lft + h/2.0, rt, f, &err2, p3, p4, p5);
      *err = err1 + err2;             /* sum the real errors, */
      return (int5);                  /* and return the integral. */
   }
}

/* Adaptive quadrature interface.
   Called with:
        lft:          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 p1, p2, p3;                 /* Initial function evaluations. */
   double h;                          /* Width of this interval. */

   h = rt - lft;                      /* Get the width of the interval. */

   p1 = (*f)(lft);                    /* Pre-compute the initial function */
   p2 = (*f)(lft + h/2.0);            /* values. */
   p3 = (*f)(rt);

                              /* Call adaptive2 to do the real work. */
   return (adaptive2 (lft, rt, f, err, p1, p2, p3));
}

/* End of File */