An automatic integrator using Trapezoidal rule

How would you know how many segments to use in a Trapezoidal rule of integration to get an accurate value of the integral?  This can be done by applying the Trapezoidal rule for 1 segment rule, then 2 segment rule, followed by 4 segment rule and so on.  As soon as the absolute relative approximate error (page 5-6) between the consecutive answers becomes less than the pre-specified tolerance chosen by the user, you would have your integral within the accuracy you desired.

Here is a MATLAB program that does that for you.  The MATLAB program that can be downloaded at http://nm.mathforcollege.com/blog/trapezoidal_rule_automatic.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://nm.mathforcollege.com/blog/html/trapezoidal_rule_automatic.html

% Simulation : Using Trapezoidal rule as an automatic integrator

% Language : Matlab 2007a

% Authors : Autar Kaw, http://nm.mathforcollege.com

% Mfile available at
% http://nm.mathforcollege.com/blog/trapezoidal_rule_automatic.m

% Last Revised : October 12, 2008

% Abstract: This program uses multiple-segment Trapezoidal
% rule to integrate f(x) from x=a to x=b within a pre-specified tolerance

clc
clear all

disp(‘This program uses multiple-segment Trapezoidal rule as an automatic integrator’)
disp(‘to integrate f(x) from x=a to x=b within a pre-specified tolerance’)
disp(‘ ‘)
disp(‘Author: Autar K Kaw.’)
disp(‘http://autarkaw.wordpress.com’)
disp(‘http://nm.mathforcollege.com’)
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=20000;
tolerance=0.005;
f=inline(‘exp(x)*sin(2*x)’);

% SIMULATION
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(‘  ‘)

% 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=2;
while (ea>tolerance) & (n<=nmax)
h=(b-a)/n;
% Keeping track of used number of segments
nused=n;
current_integral=0;
for i=1:1:n-1
current_integral=current_integral+f(a+i*h);
end
current_integral=2*current_integral+f(a)+f(b);
current_integral=(b-a)/(2*n)*current_integral;
% 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

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)
if (ea>tolerance)
disp(‘  ‘)
disp(‘  ‘)
disp(‘     NOTE: The value of integral is not within the pre-specified tolerance’)
end
disp(‘  ‘)

This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://nm.mathforcollege.com.

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.