Clear Filters
Clear Filters

Why is my Simpson's 3/8 code not producing the correct values?

11 views (last 30 days)
Hi,
I am working on some simple numerical quadrature methods, and I am trying to get the Simpson's 3/8 method to work using the following code:
function I = SimpThreeEight(f, a, b, n)
h = (b-a)/n;
S =feval(f,a);
for i = 1 : 3: n-1
x(i) = a + h*i;
S = S + 3*feval(f, x(i));
end
for i = 2 : 3: n-2
x(i) = a + h*i;
S = S + 3*feval(f, x(i));
end
for i = 3 : 3: n-3
x(i) = a + h*i;
S = S + 2*feval(f, x(i));
end
S = S + feval(f, b);
I = 3*h*S/8;
However, using test data with know values my approximations are all off. Can someone identify why?
Thanks!
  1 Comment
Zhengyu Pan
Zhengyu Pan on 24 Dec 2014
use n-1 for all three of them( see my code below, mine is a little bit different I use a: 3*h : b-h instead of 1:3: n-1, but same thing.

Sign in to comment.

Accepted Answer

Roger Stafford
Roger Stafford on 23 Nov 2013
The upper limits in the first two for-loops are not right. They should be:
for i = 1:3:n-2
x(i) = a + h*i;
S = S + 3*feval(f, x(i));
end
for i = 2:3:n-1
x(i) = a + h*i;
S = S + 3*feval(f, x(i));
end
Your second for-loop: "for i = 2:3:n-2" incorrectly leaves out the term with 3*feval(f,x(n-1)).

More Answers (2)

bym
bym on 23 Nov 2013
your first loop should start at 2, second at 3, third at 4. Set end points for all loops at n-1. Make sure n is divisible by 3 (i.e. 6,9,12 etc)

Zhengyu Pan
Zhengyu Pan on 24 Dec 2014
Edited: Zhengyu Pan on 24 Dec 2014
This is from my final, so code is way long than it should be :D
function compSimpson(f,a,b)
%f is the integrand function
% a,b are endpoints of the interval
% c is true value of the interval from a to b, calculate by hand
c =quadgk(f,a,b); % use c to calculate the integral value
int_approx=zeros; % make the variables we wanted as zero matrixs
hMatrix=zeros;
error=zeros;
s=zeros;
m=1:10; % number of panels
disp('integal area of f =')
disp(c)
disp('-------------------------------------------') % my fancy chart
disp([' m h ', ' approx', ' error slope']) % not practical,
disp('-------------------------------------------') % but fancy
for n=1:1:10
h = (b-a)/n; % length of h, the panel
sum = 0; % initial of sum 2nd, 5th 8th… term without coeffcient
sumtwo = 0; % initial of sum 3rd,6th,…. term without coeffcient
sumthree= 0; % initial of sum 4rd,7th,…. term without coeffcient
for k = a+h: 3*h:b-h
sum = sum + f(k); % sum of "y(2i-1)
end% end for k
for j= a+2*h:3*h:b-h
sumtwo = sumtwo + f(j); % sum of "y(2i)"
end % for j
for w= a+3*h: 3*h: b-h
sumthree=sumthree+f(w);
end
int_approx(n, 1) = (3*h/8)*(f(a)+f(b)+3*sum+3*sumtwo+ 2*sumthree); % approx of interval
% using Composite
% Simpson's Rule of 3/8
hMatrix(n,1)= h; % make a all h as
% matrix
error(n,1)=abs((3*h/8)*(f(a)+f(b)+3*sum+3*sumtwo+ 2*sumthree)-c); % all error terms
end % end for n
ssum=0;
for n= 2:10
s(1,1)=0;
s(n,1)= log(error(n,1)/ error(n-1,1) ) /log(hMatrix(n,1)/hMatrix(n-1,1));
% matrix for slop
ssum= ssum+ s(n,1); % sum of the the slopes
end % end for n again
averageOfSlope= log(error(9,1)/error(1,1))/log(hMatrix(9,1)/hMatrix(1,1)) ;
figure
loglog(hMatrix, error,'>-') % standand figure
%axis tight % just make the graph look good
xlabel('h') % x-axis label
ylabel('error') % y-axis label
T=evalc('f'); % transfer function to a readable string
legend(T,'Location','northwest') % to let us know what graph is that
legend('boxoff') % practical speaking, we don't need the box
%p = polyfit(hMatrix,error,1);
m =m';
%C=evalc('s(9,1)');
disp([m,hMatrix, int_approx, error, s]) % to show my fancy chart
disp('-------------------------------------------')
disp(' ')
disp(['the average of the slope goes to ',num2str(averageOfSlope)])
end % end for the function
EDU>> f=@(x) sin(pi*x)
f =
@(x)sin(pi*x)
EDU>> compSimpson(f,0,1) integal area of f = 0.6366
----------------------------------------------------------
m h approx error slope
----------------------------------------------------------
1.0000 1.0000 0.0000 0.6366 0
2.0000 0.5000 0.5625 0.0741 3.1025
3.0000 0.3333 0.6495 0.0129 4.3124
4.0000 0.2500 0.6127 0.0239 -2.1457
5.0000 0.2000 0.6211 0.0155 1.9518
6.0000 0.1667 0.6373 0.0006 17.4724
7.0000 0.1429 0.6287 0.0080 -16.3520
8.0000 0.1250 0.6305 0.0061 1.9866
9.0000 0.1111 0.6367 0.0001 33.2404
10.0000 0.1000 0.6327 0.0039 -32.9430
-------------------------------------------
the average of the slope goes to 3.897
comment:
for this problem, we can conclude that simpson’s 3/8 rule is accurate if m(number of panels) is multiple of 3. Moreover, the order of the error formula for composite Simpson’s 3/8 rule is about 4

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!