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

% Last Revised : October 8, 2008

% Abstract: This program shows use of multiple segment Trapezoidal
% rule to integrate 1/sqrt(x) from x=0 to b, b>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.
% b  = Upper limit of integration
% m = Maximum number of segments is 2^m
b=9;
m=14;


% SIMULATION
fprintf('\nFinding the integral of 1/sqrt(x) with limits of integration as x=0 to x=%g',b)

% EXACT VALUE OF INTEGRAL
% integrand 1/sqrt(x)
syms x
f=1/sqrt(x);
a=0;
valexact=double(int(f,x,a,b));
fprintf('\n\nExact value of integral = %f',valexact)
disp( '  ')

f=inline('1/sqrt(x)');
%finding value of the integral using 16,...2^m segments
for k=4:1:m
    n=2^k;
    h=(b-a)/n;
    sum=0;
        for i=1:1:n-1
            sum=sum+f(a+i*h);
        end
% See below how f(a) is not added as f(a)=infinity.  Instead we 
% use a value of f(a)=0.  How can we do that? Because as per integral calculus, 
% using a different value of the function at one point or 
% at finite number of points does not change the value of the 
% integral.
    sum=2*sum+0+f(b);
    sum=(b-a)/(2*n)*sum;
    fprintf('\nApproximate value of integral =%f with %g segments',sum,n)    
end 
disp('  ')
