MATLAB Answers

0

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

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

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.

로그인 to comment.

답변 수: 3

Roger Stafford 님의 답변 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)).

  댓글 수: 1

Thanks! I figured it was something simple.

로그인 to comment.


bym 님의 답변 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)

  댓글 수: 0

로그인 to comment.


Zhengyu Pan 님의 답변 24 Dec 2014
Zhengyu Pan 님이 편집함. 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

  댓글 수: 0

로그인 to comment.



Translated by