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

조회 수: 9(최근 30일)
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
Zhengyu Pan
Zhengyu Pan 2014년 12월 24일
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.

댓글을 달려면 로그인하십시오.

채택된 답변

Roger Stafford
Roger Stafford 2013년 11월 23일
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
Jerome
Jerome 2013년 11월 23일
Thanks! I figured it was something simple.

댓글을 달려면 로그인하십시오.

추가 답변(2개)

bym
bym 2013년 11월 23일
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 2014년 12월 24일
편집: Zhengyu Pan 2014년 12월 24일
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

범주

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!

Translated by