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