image thumbnail
from Comparison analysis of numerical intergration methods by Sulaymon Eshkabilov
Comparison analysis of numerical integration methods, viz.trapezoid, Simpsons Rule, mid-point, etc.

NUM_quadrature_compare.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                Written by Sulaymon L. ESHKABILOV, Ph.D
%                         October, 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%% NUM_quad_COMPARE.m
%{
 Comparison analysis of numerical intergration methods, viz.    trapezoid, Composite trapezoid, Simpson's Rule, Composite Simpson's rule, Mid-point Rule, Composite Mid-point Rule, etc.
%}
clear all;
clc;
% Example # 1. Polynomial integration
disp('%--------------------------------------%')
display('Example # 1. f=3*t^2-4*t+2')
disp('%--------------------------------------%')
f1=inline('3*t.^2-4*t+2');
a=0; b=1;
% Analytical solution of the given problem:
syms t1 t2
FF=3*t1^2-4*t1+2; FFi=int(FF); t1=0; % Lower limit calc's
FF1=subs(FFi);
FF=3*t2^2-4*t2+2; FFi=int(FF); t2=1; % Upper limit calc's 
FF2=subs(FFi);
FFsol=FF2-FF1; 
fprintf('ANALYTIC Solution"s numerical value is: %5.5f \n',... 
    FFsol)
%% Trapezoid rule: Calls a function file called trapezoid.m 
% For help: type in help trapezoid or doc trapezoid
disp('%--------------------------------------%')
display('Trapezoid method:')
F_trap=trapezoid(f1, b, a);
 
% ABSOLUTE and RELATIVE percentage errors of the trapezoid rule:
Trap_errABS=abs(FFsol-F_trap);  
Trap_errREL=abs(FFsol-F_trap)/(abs(FFsol))*100; 
fprintf('Trapezoid method"s Abs. error: %5.5f \n', Trap_errABS)
fprintf('Trapezoid method"s Rel. error: %5.5f \n', Trap_errREL)

%% COMPOSITE TRAPEZOID rule. Calls composite_trapezoid.m.
% For help: type in help composite_trapezoid or doc 
% composite_trapezoid in the command window
disp('%--------------------------------------%')
disp('Composite trapezoid method: ')
F_comTRAP=composite_trapezoid(f1,  b, a, 16);
% ABSOLUTE & RELATIVE percentage errors of composite trapezoid: 
comTRAP_errABS=abs(FFsol-F_comTRAP(end)); 
comTRAP_errREL=abs(FFsol-F_comTRAP(end))/(abs(FFsol))*100; 
fprintf('Composite trapezoid method"s Abs. error:  %5.5f\n', ... 
    comTRAP_errABS)
fprintf('Composite trapezoid method"s Rel. error:  %5.5f\n', ... 
comTRAP_errREL)
 
%% SIMPSON's RULE. Calls a function file called simpsons_rule.m. 
% For help about the m-file enter >> help simpsons_rule or 
% doc simpsons_rule in the command window
disp('%--------------------------------------%')
display('SIMPSONs RULE method: ')
F_SIM=simpsons_rule(f1, b, a);
% ABSOLUTE and RELATIVE percentage errors of Simpson's method
 
SIM_errABS=abs(FFsol-F_SIM); disp(SIM_errABS); 
SIM_errREL=abs(FFsol-F_SIM)/(abs(FFsol))*100; 
fprintf('Simpsons method"s Abs. error:  %5.5f \n', SIM_errABS)
fprintf('Simpsons method"s Rel. error:  %5.5f \n', SIM_errREL)
 
%% COMPOSITE SIMPSON's method. Calls COMP_simpsons_rule.m
% For HELP, type in >> help COMP_simpsons_rule.m or doc 
%       doc COMP_simpsons_rule.m
disp('%--------------------------------------%')
display('Composite Simpmson"s Rule')
F_comSIM=COMP_simpsons_rule(f1,  b,a, 16);
% ABSOLUTE and RELATIVE percentage errors of the composite simpson's method
comSIM_errABS=abs(FFsol-F_comSIM); 
comSIM_errREL=abs(FFsol-F_comSIM)/(abs(FFsol))*100; 
fprintf('Composite Simpsons method"s Abs. error: %5.5f \n', ... 
    comSIM_errABS)
fprintf('Composite Simpsons method"s rel. error: %5.5f \n', ... 
comSIM_errREL)
 
%% MIDPOINT RULE Method calls midpoint_rule.m
disp('%--------------------------------------%')
disp('Midpoint rule method')
F_midP=midpoint_rule(f1, b, a);
% ABSOLUTE and RELATIVE percentage errors of the midpoint rule: 
midP_errABS=abs(FFsol-F_midP);  
midP_errREL=abs(FFsol-F_midP)/(abs(FFsol))*100; 
fprintf('Midpoint rule method"s Abs. error: %5.5f \n', ...
    midP_errABS)
fprintf('Midpoint rule method"s Rel. error: %5.5f \n', ... 
midP_errREL)
 
%% COMPOSITE MIDPOINT Rule calls comp_midpoint.m
disp('%--------------------------------------%')
display('Composite MIDpoint method')
F_comMP=comp_midpoint_rule(f1,  b, a,16);
% ABS. and REL. percentage errors of composite midpoint method
comMP_errABS=abs(FFsol-F_comMP);  
comMP_errREL=abs(FFsol-F_comMP)/(abs(FFsol))*100; 
fprintf('Composite midpoint method"s Abs. error: %5.5f \n', ... 
    comMP_errABS)
fprintf('Composite midpoint method"s Abs. error: %5.5f \n', ... 
comMP_errREL)
 
%%  QUAD MATLAB function uses adaptive Simpson quadrature.
% Q = QUAD(FUN,A,B) approximates the integral of scalar-valued
% function FUN from A to B to within an error of 1.e-6 using     % recursive adaptive Simpson quadrature. FUN is a function       % handle. The function Y=FUN(X) should accept a vector argument % X and return a vector result Y, the integrand evaluated at    % each element of X.
disp('%---------------------------------------------------%')
display('QUAD MATLAB built function tool:');
F_MQ=quad(f1, a, b); 
fprintf('Final value of Integration is: %25.5f\n',F_MQ);
% NB: watch-out!!! position of integral limits a & b for quad 
 
% ABSOLUTE and RELATIVE percentage errors of the MATLAB QUAD 
MQ_errABS=abs(FFsol-F_MQ); 
MQ_errREL=abs(FFsol-F_MQ)/(abs(FFsol))*100; 
fprintf('MATLAB Quad approx. method"s Abs. error: %5.5f \n', ... 
    MQ_errABS)
fprintf('MATLAB Quad approx. method"s Rel. error: %5.5f \n', ... 
MQ_errREL)
 
 
%% Example # 2
clear all
disp('%--------------------------------------%')
disp('Example # 2. f=2*sin(2*t)+tan(t) ')
disp('%--------------------------------------%')
f2=@(t)(2*sin(2*t)+tan(t));
a=30; b=45;
 
% Analytic solution of the given problem
syms t1 t2
FF=2*sin(2*t1)+tan(t1); FFi=int(FF); t1=deg2rad(30); FF1=subs(FFi);
FF=2*sin(2*t2)+tan(t2); FFi=int(FF); t2=deg2rad(45); FF2=subs(FFi);
disp('%--------------------------------------%')
display('ANALYTIC Solution"s numerical value')
FFsol=FF2-FF1; 
fprintf('ANALYTIC Solution"s numerical value: %25.5f\n', FFsol)
% trapezoid
disp('%--------------------------------------%')
tic;
display('Output of the trapezoid method  = ')
F_trap=trapezoid(f2, b, a, 'trig');
Time_trapezoid=toc;
 
TRAP_errABS=abs(FFsol-F_trap); 
TRAP_errREL=abs(FFsol-F_trap)/(abs(FFsol))*100; 
fprintf('trapezoid method"s Abs. error: %5.5f \n', TRAP_errABS)
fprintf('trapezoid method"s Rel. error: %5.5f \n', TRAP_errREL)
 
% COMPOSITE trapezoid.
disp('%--------------------------------------%')
tic;
disp('OUT_composite_trapezoid = ')
F_CTRAP=composite_trapezoid(f2, a, b, 56, 'trigonometric function');
Time_composite_trap=toc;
 
CTRAP_errABS=abs(FFsol-F_CTRAP(end)); 
CTRAP_errREL=abs(FFsol-F_CTRAP(end))/(abs(FFsol))*100; 
fprintf('Composite trapezoid method"s Abs. error: %5.5f \n', ... 
    CTRAP_errABS)
fprintf('Composite trapezoid method"s Rel. error: %5.5f \n', ... 
CTRAP_errREL)
 
% SIMPSON's RULE.
disp('%--------------------------------------%')
tic;
disp('SIMPSONs Rule')
F_SR=simpsons_rule(f2, a,b, 'trig. function');
Time_simpsonsRULE=toc;
 
SR_errABS=abs(FFsol-F_SR); 
SR_errREL=abs(FFsol-F_SR)/(abs(FFsol))*100; 
fprintf('Simpson"s rule method"s Abs. error:  %5.5f \n', ...
    SR_errABS)
fprintf('Simpson"s rule method"s Rel. error:  %5.5f \n', ... 
SR_errREL)
 
 
% COMPOSITE SIMPSON's Rule. 
% NB: The number of steps in this method has to be even.
disp('%--------------------------------------%')
tic;
display('COMPOSITE SIMPSON"s Rule');
F_CSR=COMP_simpsons_rule(f2, a, b,56, 'trig. function' );
Time_comp_simpsonsRULE=toc;
 
Csr_errABS=abs(FFsol-F_CSR);  
Csr_errREL=abs(FFsol-F_CSR)/(abs(FFsol))*100; 
fprintf('Composite Simpson"s rule Abs. error: %5.5f\n', ... 
    Csr_errABS)
fprintf('Composite Simpson"s rule Rel. error: %5.5f\n', ... 
Csr_errREL)
 
% MIDPOINT Rule method. 
disp('%--------------------------------------%')
tic
display('MIDPOINT Rule Method');
F_MP=midpoint_rule(f2, a, b,56, 'trig. function' );
Midpoint_RULE=toc;
 
MP_errABS=abs(FFsol-F_MP); 
MP_errREL=abs(FFsol-F_MP)/(abs(FFsol))*100; 
fprintf('Midpoint rule method"s Abs. error: %5.5f \n', ...
    MP_errABS)
fprintf('Midpoint rule method"s Rel. error: %5.5f \n', ...
MP_errREL)
 
% COMPOSITE MIDPOINT method. 
disp('%--------------------------------------%')
tic
display('Composite Midpoint rule METHOD ');
F_Cmp=comp_midpoint_rule(f2, a,b,56, 'trig. function ');
COMP_midpoint_rule=toc;
 
Cmp_errABS=abs(FFsol-F_Cmp);  
Cmp_errREL=abs(FFsol-F_Cmp)/(abs(FFsol))*100; 
fprintf('Composite midpoint rule"s Abs. error: %5.5f\n', ...
Cmp_errABS)
fprintf('Composite midpoint rule"s Rel. error: %5.5f\n', ...
Cmp_errREL)
 
%%  MATLAB QUAD:   
disp('%--------------------------------------%')
tic
display('QUAD MATLAB builtin numeric eval of integrals');
F_MQ=quad(f2, deg2rad(a), deg2rad(b)); disp(F_MQ);
Time_QUAD=toc;
MQ_errABS=abs(FFsol-F_MQ); 
MQ_errREL=abs(FFsol-F_MQ)/(abs(FFsol))*100; 
fprintf('QUAD MATLAB function"s Abs. error: %5.5f\n ', ...
    MQ_errABS)
fprintf('QUAD MATLAB function"s Rel. error: %5.5f\n ', ... 
MQ_errREL)

disp('Computation time of quadrature rules:')
fprintf('MATLAB quad"s time for ex2 is: %5.5f\n', Time_QUAD)
fprintf('Trapezoid time for ex2 is: %5.5f\n', Time_trapezoid)
fprintf('Comp. Trapezoid"s time for ex2 is: %5.5f\n', ... 
    Time_composite_trap)
fprintf('Simpson"s Rule time for ex2 is: %5.5f\n', ... 
Time_simpsonsRULE)
fprintf('Composite Simpson"s time for ex2 is: %5.5f\n', ... 
Time_comp_simpsonsRULE)
 

Contact us