#!/usr/bin/env python

from math import *
from adaptive_simpson import adaptive_simpson

#-----------------------------------------------------------------------------
# Jonathan Senning <jonathan.senning@gordon.edu>
# February 20, 2007
#
# $Id$
#
# This program verifies that the adaptive Simpson's rule scheme described
# in "Numerical Mathematics and Computing" (5th edition) by Cheney and
# Kincaid, Brooks/Cole, 2004.
#
# The key assumption in the algorithm is that the 4th derivative of the
# integrand f(x) has the same value everywhere on an interval.  This is not
# verified by the algorithm, but presumably as the interval size is reduced
# it becomes more reasonable.
#
# This program computes the integral of f(x) on [a,b] using the adaptive
# Simpson's rule and computes the difference between the absolute value of the
# actual error and the specified tolerance. If the assumption above is correct
# then the specified tolerance should be satisifed.

#-----------------------------------------------------------------------------

# Define the interval of integration

a =  0.0
b = 10.0

# Define the integrand and its antiderivative (used to compute the
# exact answer.

def f( x ): return exp( -x ) * sin( x )
def exact( x ): return 0.5 * ( 1 - exp( -x ) * ( sin( x ) + cos( x ) ) )

#-----------------------------------------------------------------------------

# Evaluate the integral with successively smaller tolerance values and report
# both the tolerance and the "tolerance error."  If this last quantity is
# negative, it indicates that the actual error is greater than the specified
# tolerance.

print "  Integral       Error     Tolerance    Message"
print "------------ ------------ ------------ -------------------------------"

tol = 1.0

for i in xrange( 30 ):
    int = adaptive_simpson( f, a, b, tol )
    err = int - exact( b )
    if tol - abs( err ) < 0:
        message = 'Error exceeds tolerance!'
    else:
        message = ''
    print "%12.5e %12.5e %12.5e %s" % ( int, err, tol, message )
    tol = tol / 2.0

#-----------------------------------------------------------------------------

# End of File