#!/usr/bin/octave -q

%-----------------------------------------------------------------------------
% 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.

function y = f(x)
  y = exp(-x) * sin(x);
endfunction

function y = exact(x)
  y = 0.5 * ( 1 - exp(-x) * ( sin(x) + cos(x) ) );
endfunction

%-----------------------------------------------------------------------------

% 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.

fprintf( '  Tolerance  ToleranceErr\n' );
fprintf( '------------ ------------\n' );

tol = 1.0;

for i = 1 : 30
    err = adaptSimp( @f, a, b, tol ) - exact(b);
    fprintf( '%12.5e %12.5e\n', tol, tol - abs(err) );
    tol = tol / 2.0;
endfor

%-----------------------------------------------------------------------------

% End of File