Simpson's 3/8 Rule
15 views (last 30 days)
Show older comments
Could someone please take a look at my code. When I execute this code it says: ??? Index exceeds matrix dimensions. Plus it gave me the wrong value. I don't know what went wrong. Thank you in advance.
% Inputs
% This is the only place in the program where the user makes the changes
% based on their wishes
% f(x), the function to integrate
f= @(x) 0.2 + 25*x - 200*x.^2 + 675*x.^3 - 900*x.^4 + 400*x.^5 ;
% a, the lower limit of integration
a= 0 ;
% b, the upper limit of integration
b= 0.8 ;
% n, the maximum number of segments
n=5 ;
%**********************************************************************
% Displays what inputs are being used
disp(sprintf('\n\n********************************Input Data**********************************\n'))
disp(sprintf(' f(x), integrand function'))
disp(sprintf(' a = %g, lower limit of integration ',a))
disp(sprintf(' b = %g, upper limit of integration ',b))
disp(sprintf(' n = %g, number of subdivision',n))
format short g
% Determine the exact solution
exact = quad(f,a,b) ;
% Determine number of different segments to consider
nstep = floor(log2(n)) ;
for i=0:nstep-1
% Determine number of segments and h
NN(i+1)=2^(i+1) ;
h=(b-a)/NN(i+1) ;
% Use the rule to find the approximation
integral = f(a) + f(b) ;
for j=1:2:NN(i+1)-1
integral = integral + 3*f(a+h*j) ;
end
for j=2:2:NN(i+1)-2
integral = integral + 3*f(a+h*j) ;
end
for j=3:2:NN(i+1)-3
integral = integral + 2*f(a+h*j) ;
end
integral = integral * 3 * h/8 ;
YY(i+1)=integral ;
% Compute Errors
Et(i+1)=exact-integral ;
Etabs(i+1)=abs((integral-exact)/exact) ;
if(i > 0)
Ea(i+1)=YY(i+1)-YY(i) ;
Eaabs(i+1)=abs((YY(i+1)-YY(i))/YY(i)) ;
SD(i+1)=floor((2-log10(Eaabs(i+1)/0.5))) ;
if(SD(i+1)<0)
SD(i+1)=0
end
else
Ea(1)=0 ;
Eaabs(1)=0 ;
SD(1)=0 ;
end
end
% The following builds and displays a table of values
disp(sprintf('\n\n************************Table of Values******************************\n'))
disp(' Approx True Relative Approx Rel Appr Sig ')
disp(' n Integral Error True Error Error Error Digits ')
disp('-----------------------------------------------------------------------')
for i=1:nstep
if(i > 1)
if(exact || YY(i) > 0)
disp(sprintf('%4i %+1.3e %+1.3e %+1.3e %+1.3e %+1.3e %2i',NN(i),YY(i),Et(i),Etabs(i),Ea(i),Eaabs(i),SD(i) ))
else
disp(sprintf('%4i %+1.3e %+1.3e n/a %+1.3e n/a n/a',NN(i),YY(i),Etabs(i),Ea(i)))
end
else
disp(sprintf('%4i %+1.3e %+1.3e %+1.3e n/a n/a n/a',NN(i),YY(i),Et(i),Etabs(i)))
end
end
disp('-----------------------------------------------------------------------')
% This function detects information about your
% screensize and tries to then place/size the graphs accordingly.
scnsize = get(0,'ScreenSize');
% Graph: The following code sets up the graphical depiction of the method.
x(1)=a ;
y(1)=f(a) ;
hold on
for i=1:n
x(i+1)=a+i*h ;
y(i+1)=f(x(i+1)) ;
end
for i=1:2:n
bg = i ;
ed = i + 2 ;
coeffs = polyfit(x(bg:ed), y(bg:ed), 2);
poly_x1 = [x(bg):(x(ed) - x(bg))/1000:x(ed)];
poly_y1 = coeffs(1)*poly_x1.^2 + coeffs(2)*poly_x1 + coeffs(3);
poly_x1 = [poly_x1(1) poly_x1 poly_x1(end)];
poly_y1 = [0 poly_y1 0];
fill(poly_x1, poly_y1, 'y')
end
xrange=a:(b-a)/1000:b;
plot(xrange,f(xrange),'k','Linewidth',2)
title('Integrand function and Graphical Depiction of Simpson''s 3/8 Rule')
0 Comments
Accepted Answer
Oleg Komarov
on 23 Mar 2012
The problem is in the last LOOP, you set at the last iteration:
ed = i + 2;
When i = 5, ed = 7 but x has only 6 elements.
EDIT Also note that on FEX just recently somebody submitted: http://www.mathworks.com/matlabcentral/fileexchange/33493-simpsons-13-and-38-rules
3 Comments
Oleg Komarov
on 24 Mar 2012
Try:
for i=2:2:n-1
bg = i-1 ;
ed = i + 1 ;
coeffs = polyfit(x(bg:ed), y(bg:ed), 2);
...
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!