Integral of a Real-valued Function over an Interval.
                                       
DESCRIPTION:

   Approximates the integral of a real-valued function over a given
   interval, and estimates the absolute error in the approximation.
   
USAGE:

integrate(f, lower, upper, subdivisions=100,
     rel.tol=.Machine$double.eps^.25,
     abs.tol=rel.tol, keep.xy=F, aux = NULL, ...)

REQUIRED ARGUMENTS:

   f
          a real-valued S-PLUS function of the form `f(x,)', where x is
          the variable of integration. Users can accumulate information
          through attributes of the function value. If the attributes
          include any additional arguments of f, the next call to f will
          use the new values of those arguments. f should be vectorized:
          the output of f(x,...) should be the length of x and
          f(x,...)[i] should be the same as f(x[i],...).
          
   lower, upper
          the lower and upper limits of integration. Infinite limits are
          allowed.
          
OPTIONAL ARGUMENTS:

   subdivisions
          an upper bound on the number of subdivisions of the original
          interval.
          
   rel.tol
          relative error tolerance. The method attempts to compute an
          approximation Q to the true integral I such that `abs(I-Q) <=
          max(abs.tol,rel.tol*abs(Q)'. The algorithm will not proceed if
          `abs.tol <= 0' and `rel.tol <
          max(50*.Machine$double.eps,.5e-28)'.
          
   abs.tol
          absolute error tolerance. The method attempts to compute an
          approximation Q to the true integral I such that `abs(I-Q) <=
          max(abs.tol,rel.tol*abs(Q))'.
          
   keep.xy
          tells whether or not to accumulate the points used by integrate
          and the corresponding function values as x and y components in
          the output.
          
   aux
          restart information from a previous call with smaller max.sub
          and / or larger tolerance. If aux is not NULL, then the
          algorithm is restarted using the information in aux.
          
   ...
          additional arguments for f.
          
VALUE:

   a list containing the following elements :
          
   integral
          the approximate integral of f from lower to upper.
          
   abs.error
          an upper bound on the magnitude of the absolute error in the
          approximation.
          
   f.last
          if one or more singularities are encountered, this is a matrix
          whose first row lists the last points at which the f was
          evaluated and whose second row lists the corresponding values
          of f.
          
   subdivisions
          the number of subintervals to which the quadrature rule has
          been applied.
          
   message
          a statement of the reason for termination.
          
   aux
          a list that can be used for restarting the algorithm allowing
          for an increase in the number of subintervals and/or a
          decreased tolerance. Includes the attributes returned after the
          last function call, as well as the function and limits of
          integration.
          
   call
          a copy of the call to integrate.
          
   x
          the vector of points at which f is evaluated by integrate (in
          order). Present only if keep.xy = T.
          
   y
          the function values associated with x. Present only if keep.xy
          = T.
          
NOTE:

   Function values should be computed in C or Fortran within the outer
   S-PLUS function for greater efficiency. When restarting, it is not
   necessary to supply the function, limits of integration, or
   appropriate attribute values since they are saved in aux.
   
METHOD:

   Implements adaptive 15-point Gauss-Kronrod quadrature based on the
   Fortran function dqage and dqagie from QUADPACK (Piessens et al. 1983)
   in NETLIB (Dongarra and Grosse 1987).
   
REFERENCES:

   Piessens, R., deDoncker-Kapenga, E., Uberhuber, C., and Kahaner, D.
   (1983). QUADPACK: A Subroutine Package for Automatic Integration.
   Springer, Berlin. Dongarra, J. J., and Grosse, E. (1987). Distribution
   of mathematical software via electronic mail, Communications of the
   ACM 30 , pp. 403-407.
   
EXAMPLES:

# Example 1:
# erf gives an approximate value for the error function
erf <- function(b, tol = .Machine$double.eps^0.5)
   {
     integrand <- function( x) {(2/sqrt(pi))*exp(-x^2)}
     integrate(integrand, lower = 0, upper = b, tol.abs = tol)$integral
   }
erf(1)
# Example 2:
# this application approximates the value of pi
integrand <- function(x) {1/{(x+1)*sqrt(x)}}
integrate(integrand, lower = 0, upper = Inf)
# Example 3:  A double integral demo to evaluate the
# double integral: integral(0,5) integral(0,t) f(x) dx dt
# where f(x) is x^2.  Result should be 5^4/12, or 52.08333.
double.integral.demo.fn <-
function()
{
        fun <- function(upper, integrand)
        {
                unlist(lapply(upper, function(upper, integrand)
                {
                        integrate(f = integrand,
                                  lower = 0,
                                  upper = upper)$integral
                }
                , integrand))
        }
        return(integrate(f = fun,
                         lower = 0,
                         upper = 5,
                         integrand = function(x)
                                        {
                                                x^2
                                        }
                )$integral)
}