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

조회 수: 9(최근 30일)
Jerome 2013년 11월 22일
댓글: Zhengyu Pan 2014년 12월 24일
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 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 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 2013년 11월 23일
Thanks! I figured it was something simple.

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

### 추가 답변(2개)

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)
##### 댓글 수: 0표시숨기기 이전 댓글 수: -1

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

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!