simpdemo.m

% Simpson's rule integration demo

fprintf('\n Simpsons 1/3 rule integration demo\n');

n = 2;
olderror = 1;
fprintf('\n  n        integral       error      fract prev error    flops\n\n');
while n < 32768
  flops(0);         % floating point count to zero
  t = clock;        % current time

  a = 0;
  b = pi;
  h = (b-a)/n;
  x = a:h:b;
  y = sin(x);
  simpint = simp0(y,h);

  error = 2.0 - simpint;  % exact result of integral is 2.0
  ratio = error/olderror;
  olderror = error;
  fprintf('%6.0f %14.9f %14.5e %10.3f % 14.0f\n',n,simpint,error,ratio,flops);

  n = n*2;
end;

fprintf('\nelapsed time: %4.2f\n\n',etime(clock,t));