% 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('  ')

This program shows the convergence of getting the value of 
an improper integral using multiple segment Trapezoidal rule
Author: Autar K Kaw.  autarkaw.wordpress.com

Finding the integral of 1/sqrt(x) with limits of integration as x=0 to x=9

Exact value of integral = 6.000000  

Approximate value of integral =4.904246 with 16 segments
Approximate value of integral =5.225408 with 32 segments
Approximate value of integral =5.452337 with 64 segments
Approximate value of integral =5.612757 with 128 segments
Approximate value of integral =5.726182 with 256 segments
Approximate value of integral =5.806382 with 512 segments
Approximate value of integral =5.863092 with 1024 segments
Approximate value of integral =5.903191 with 2048 segments
Approximate value of integral =5.931546 with 4096 segments
Approximate value of integral =5.951596 with 8192 segments
Approximate value of integral =5.965773 with 16384 segments  

Published with MATLAB® 7.3