%% decrease_h
%
%  Discussion:
%
%    Look at how trapezoid quadrature error behaves as we decrease h.
%
%    Integrate y=exp(x) from -2^k to +2^k, for k = 0, -1, ..., -10.
%    Expect error approximately c*h^3.
%    Hence error should decrease by factor 1/8 each step.
%
format long

fprintf ( 1, '\n' );
fprintf ( 1, 'decrease_h:\n' );
fprintf ( 1, '  Look at how trapezoid quadrature error behaves as we\n' );
fprintf ( 1, '  decrease interval width h\n' );

fprintf ( 1, '\n' );
fprintf ( 1, ' k     xleft     xright     int         quad       error     rate\n' );
fprintf ( 1, '\n' );

e = 0.0;
twok = 1.0;
for k = 0 : 10
  xleft = - twok;
  xright = + twok;
  i = i_exact ( xleft, xright );
  q = ( xright - xleft ) * 0.5 * ( f(xleft) + f(xright) );
  eold = e;
  e = abs ( i - q );
  if ( eold ~= 0.0 )
    r = e / eold;
  else
    r = 0;
  end
  fprintf ( 1, '%2d  %8f  %8f  %10.6g  %10.6g  %10.6g  %6.4f\n', ...
    k, xleft, xright, i, q, e, r );

  twok = twok / 2.0;

end

fprintf ( 1, '\n' );
fprintf ( 1, '  Normal end of execution.\n' );

function value = f(x)
  value = exp ( x );
  return
end
function value = i_exact ( a, b )
  value = exp(b) - exp(a);
  return
end

