Another improper integral solved using trapezoidal rule

In a previous post, I showed how Trapezoidal rule can be used to solve improper integrals.  The example used in the post was an improper integral with an infinite interval of integration.

In an example in this post, we use Trapezoidal rule to solve an improper integral where the integrand becomes infinite.  The integral is .  The integrand becomes infinite at x=0.  Since x=0 would be one of the points where the integrand will be sought by the multiple-segment Trapezoidal rule, we choose the value of the integrand at x=0 to be zero (any other value would do too – a better assumption would be f(h), where h is the segment width in the multiple-segment Trapezoidal rule).

Here is a MATLAB program that shows you the exact value of the integral and then compares it with the multiple-segment Trapezoidal rule.  The convergence is slow but you can integrate improper integrals using Trapezoidal rule.

The MATLAB program that can be downloaded at http://numericalmethods.eng.usf.edu/blog/trapezoidal_improper_sqrtx.m (better to download it as single quotes from the web-post do not translate correctly with the MATLAB editor).  The html file showing the mfile and the command window output is here: http://numericalmethods.eng.usf.edu/blog/html/trapezoidal_improper_sqrtx.html

% 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 post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu.

An abridged (for low cost) book on Numerical Methods with Applications will be in print (includes problem sets, TOC, index) on December 10, 2008 and available at lulu storefront.

Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.