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

6 次查看（过去 30 天）
Jerome
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
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)).

更多回答（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

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by