SIMULATION : AUTOMATIC INTEGRATOR - DOUBLING SEGMENTS TRAPEZOIDAL RULE FORMULA

Language : Matlab 2008a Authors : Autar Kaw, http://numericalmethods.eng.usf.edu Mfile available at http://numericalmethods.eng.usf.edu/blog/trapezoidal_rule_automatic.m Last Revised : March 5, 2009 Abstract: This program uses multiple-segment Trapezoidal rule to integrate f(x) from x=a to x=b within a pre-specified tolerance

Contents

INTRODUCTION

clc
clear all

disp('This program uses special formula of multiple-segment Trapezoidal rule')
disp ('as an automatic integrator to integrate f(x) from x=a to x=b')
disp ('within a pre-specified tolerance')
disp(' ')
disp('Author: Autar K Kaw.')
disp('http://autarkaw.wordpress.com')
disp('http://numericalmethods.eng.usf.edu')
disp(' ')

This program uses special formula of multiple-segment Trapezoidal rule
as an automatic integrator to integrate f(x) from x=a to x=b
within a pre-specified tolerance
 
Author: Autar K Kaw.
http://autarkaw.wordpress.com
http://numericalmethods.eng.usf.edu
 

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=200000;
tolerance=0.000005;
f=inline('exp(x)*sin(2*x)');

DISPLAYING INPUTS

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

INPUTS
     The integrand is =exp(x)*sin(2*x)
     Lower limit of integration, a= 5.3
     Upper limit of integration, b= 10.7
     Maximum number of segments, nmax = 200000
     Pre-specified percentage tolerance, eps = 5e-006  
  

THE CODE

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

DISPLAYING OUTPUTS

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)
fprintf('\n      Time Used =%g seconds', time_used)
if (ea>tolerance)
    disp('  ')
    disp('  ')
    disp('     NOTE: The value of integral is not within the pre-specified tolerance')
end
disp('  ')

OUTPUTS
      Number of segments used  =32768
      Approximate value of integral is =19681.6
      Absolute percentage relative approximate error =1.67643e-006
      Time Used =37.624 seconds  

Published with MATLAB® 7.6