Code covered by the BSD License  

Highlights from
Composite Simpsons method for numeric integration

image thumbnail
from Composite Simpsons method for numeric integration by Sulaymon Eshkabilov
COMPOSITE SIMPSON's method for numerical calculations and analysis exercises of Numeric Integration.

F=COMP_simpsons_rule(f,lowl, upl,nsteps, varargin)
function F=COMP_simpsons_rule(f,lowl, upl,nsteps, varargin)

% HELP. COMPOSITE SIMPSON's method. Some numerical calculations
%       and analysis exercises of Numeric Integration. 
%
%       f function is given in terms of a symbolic variable x and as an
%       inline function. E.g., f=inline('x^2+2*x-2'). Also, if the function
%       f is trigonometric function, the 4th argument can be entered as
%       'trigonom' or just 'trig' or 1. X is expected to be in degrees for
%       trigonometric function evaluations. The number of steps NSTEPS has 
%       to be even.
%       upl and lowl are upper and lower limits. NB: A sequential order of 
%       limits is unnecessary to follow, 'if' conditions will take care of
%       lower and upper limits accordingly.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                Written by Sulaymon L. ESHKABILOV, Ph.D
%                         October, 2011
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%

if rem(nsteps,2)~=0
    
    display(' So re-enter new even NSTEPS!!!')
    error('NSTEPS must be even')
    
end

if nargin==4
    if upl>lowl
        hsize=(upl-lowl)/nsteps;
        x0=lowl;
        F0=f(x0);
        F1=0;F2=0;
        n=1:nsteps;
        x=x0+n.*hsize;
        Fn=f(x);
        for ii=1:nsteps/2;
            F2=F2+(hsize/3)*(4*f(x(2*ii-1)));
        end
        for jj=2:nsteps/2;
            F1=F1+(hsize/3)*(2*f(x(2*jj-2)));
        end
        F0n=(hsize/3)*(F0+Fn(end));
        F=F1+F2+F0n;
    else
        hsize=(lowl-upl)/nsteps;
        x0=upl;
        F0=f(x0);
        F1=0;F2=0;
        mm = 1:nsteps;
        x  = x0+mm.*hsize;
        Fn = f(x);
        for ii=1:nsteps/2;
            F2=F2+(hsize/3)*(4*f(x(2*ii-1)));
        end
        for  jj=2:nsteps/2;
            F1=F1+(hsize/3)*(2*f(x(2*jj-2)));
        end
        F0n=(hsize/3)*(F0+Fn(end));
        F=F1+F2+F0n;
    end
else
    if upl>lowl
        hsize=(upl-lowl)/nsteps;
        x0=lowl;
        F0=f(deg2rad(x0));
        F1=0;F2=0;
        k  = 1:nsteps;
        x  = deg2rad(x0+k.*hsize);
        Fn = f(x);
        for ii=1:nsteps/2;
            F2=F2+deg2rad(hsize/3).*(4.*f(x(2*ii-1)));
        end
        for jj=2:nsteps/2;
            F1=F1+deg2rad(hsize/3)*(2*f(x(2*jj-2)));
        end
        F0n=deg2rad(hsize/3)*(F0+Fn(end));
        F=F1+F2+F0n;
    else
        hsize=(lowl-upl)/nsteps;
        x0=upl;
        F0=f(deg2rad(x0));
        F1=0;F2=0;
        m  = 1:nsteps;
        x  = deg2rad(x0+m.*hsize);
        Fn = f(x);
        for ii=1:nsteps/2;
            F2=F2+deg2rad(hsize/3)*(4*f(x(2*ii-1)));
        end
        for jj=2:nsteps/2;
            F1=F1+deg2rad(hsize/3)*(2*f(x(2*jj-2)));
        end
        F0n=deg2rad(hsize/3)*(F0+Fn(end));
        F=F1+F2+F0n;
        
    end
end
fprintf('Final value of Integration is: %25.6f  \n', F);
return

Contact us