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