Asked by Jerome
on 22 Nov 2013

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!

Answer by Roger Stafford
on 23 Nov 2013

Accepted Answer

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)).

Answer by Zhengyu Pan
on 24 Dec 2014

Edited by 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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 1 Comment

## Zhengyu Pan (view profile)

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/107240-why-is-my-simpson-s-3-8-code-not-producing-the-correct-values#comment_257394

Sign in to comment.