% Simulation : Can I use Trapezoidal rule for an improper integral?
% Language : Matlab 2007a
% Authors : Autar Kaw, http://numericalmethods.eng.usf.edu
% mfile available at
% http://numericalmethods.eng.usf.edu/blog/trapezoidal_improper.m
% Last Revised : July 15, 2008
% Abstract: This program shows use of multiple segment Trapezoidal
% rule to integrate exp(-x) from x=a to infinity, where a>0.
clc
clear all
disp('This program shows the convergence of getting the value of ')
disp('an improper integral using multiple segment Trapezoidal rule')
disp('Author: Autar K Kaw.  autarkaw.wordpress.com')

%INPUTS.  If you want to experiment, these are the only variables  
% you should and can change.
% a = Lower limit of integration
a=1/2;

fprintf('\nFinding the integral of exp(-x) from a to infinity where a=%g',a)

% SIMULATION
syms x
% integrand exp(-x)
f=exp(-x);
valexact=double(int(f,x,a,inf));
fprintf('\n\nExact value of integral = %f',valexact)
disp( '  ')

% Transformed integrand by change of variables
ft=inline('exp(-1/t)*(-1/t^2)');
% Limits of integration
at=1/a;
b=0;

%finding value of the integral using 16, 32, 64, 128 and 256 segments
for k=4:1:8
    n=2^k;
    h=(b-at)/n;
    sum=0;
        for i=1:1:n-1
            sum=sum+ft(at+i*h);
        end
% See how f(b) is not added as f(b)=infinity.  But 
% using a different value of the function at one point or 
% finite number of points does not change the value of the 
% integral. We are substituting f(b)=0
    sum=2*sum+ft(at)+0;
    sum=(b-at)/(2*n)*sum;
    fprintf('\nApproximate value of integral =%f with %g segments',sum,n)    
end 
disp('  ')
