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!

[Zhengyu Pan]

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.

Answer 1

[Roger Stafford]

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

[Jerome]

Thanks! I figured it was something simple.

Answer 2

[bym]

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)

Answer 3

[Zhengyu Pan]

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