# trapadapt.R
# Richard Taylor 19/11/2018
# Implementation of adaptive quadrature based on the trapezoid rule.
# (A better approach would economize by re-using function values already calculated.)

# function to approximate the integral of f on [a,b] by a 2-point trapezoid rule:
trap1 <- function( f, a, b ) {
  return( 0.5 * (b-a) * (f(a) + f(b)) );
}

# ====================================================================================
# function to approximate the integral of f on [a,b] by recursively subdividing by
# trapezoids until the total error estimate is less than epsilon
trap.adapt <- function( f, a, b, epsilon ) {

 # estimate using one trapezoid:
 I1 <- trap1( f, a, b );

 # improved estimate using two trapezoids:
 I2 <- trap1( f, a, (a+b)/2 ) + trap1( f, (a+b)/2, b );

 # error estimate:
 err <- abs( I1 - I2 ) / 3;   # why 1/3? ...see class notes

 cat( "got error estimate", err, "on [", a, ",", b, "] so " );

 if (err < epsilon) {  # is error small enough?
   cat("TERMINATING RECURSION.\n");
   return( I2 );
 } else {              # if not, we subdivide further...
   cat("subdividing further...\n");
   return( trap.adapt( f, a, (a+b)/2, epsilon/2 ) + trap.adapt( f, (a+b)/2, b, epsilon/2 ) );
 }

}
# ====================================================================================

# use our method on f(x)=x^2:
# (we use an "anonymous" function here since it isn't worth giving our function a name)
I <- trap.adapt( function(x) x^2, 0, 1, epsilon=1.e-3 );

cat( sprintf( "\nThe integral evaluates to %0.12f\n", I ) );