%% SIMULATION : AUTOMATIC INTEGRATOR - DOUBLING SEGMENTS TRAPEZOIDAL RULE FORMULA
% Language : Matlab 2008a
% Authors : Autar Kaw, http://numericalmethods.eng.usf.edu
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/trapezoidal_rule_automatic.m
% Last Revised : March 5, 2009
% Abstract: This program uses multiple-segment Trapezoidal
% rule to integrate f(x) from x=a to x=b within a pre-specified tolerance

%% INTRODUCTION
clc
clear all

disp('This program uses special formula of multiple-segment Trapezoidal rule') 
disp ('as an automatic integrator to integrate f(x) from x=a to x=b') 
disp ('within a pre-specified tolerance')
disp(' ')
disp('Author: Autar K Kaw.')
disp('http://autarkaw.wordpress.com')
disp('http://numericalmethods.eng.usf.edu')
disp(' ')

%% INPUTS
% If you want to experiment, these are the only variables  
% you should and can change.
% a = Lower limit of integration
% b = Upper limit of integration
% nmax = Maximum number of segments
% tolerance = pre-specified tolerance in percentage
% f = inline function as integrand
a=5.3;
b=10.7;
nmax=200000;
tolerance=0.000005;
f=inline('exp(x)*sin(2*x)');

%% DISPLAYING INPUTS
disp('INPUTS')
func=['     The integrand is =' char(f)];
disp(func) 
fprintf('     Lower limit of integration, a= %g',a)
fprintf('\n     Upper limit of integration, b= %g',b)
fprintf('\n     Maximum number of segments, nmax = %g',nmax)
fprintf('\n     Pre-specified percentage tolerance, eps = %g',tolerance)
disp('  ')
disp('  ')

%% THE CODE
% Doing the automatic integration
% Calculating the integral using 1-segment rule
previous_integral=(b-a)/2*(f(a)+f(b));
% Initializing ea as greater than pre-specified tolerance for loop to work
ea=2*tolerance;
% Starting with 2-segments inside the while loop
n=1;
time_start=clock;
while (ea>tolerance) & (n<=nmax)
    g=(b-a)/(2*n);
    % Keeping track of used number of segments
    nused=n*2;
    additional_values=0;
        for i=1:1:n
            additional_values=additional_values+f(a+(2*i-1)*g);
        end
    current_integral=previous_integral/2+additional_values*g;
    % Calculating the absolute relative approximate error
    ea = abs((current_integral-previous_integral)/current_integral)*100;
    previous_integral=current_integral;
    % Doubling the number of segments for next estimate of the integral
    n=n*2;
end
time_end=clock;
time_used=etime(time_end,time_start);

%% DISPLAYING OUTPUTS
disp('OUTPUTS')
fprintf('      Number of segments used  =%g', nused)
fprintf('\n      Approximate value of integral is =%g',current_integral)
fprintf('\n      Absolute percentage relative approximate error =%g', ea)
fprintf('\n      Time Used =%g seconds', time_used)
if (ea>tolerance) 
    disp('  ')
    disp('  ')
    disp('     NOTE: The value of integral is not within the pre-specified tolerance')
end
disp('  ')
