%  Tests routines trap and simp for computing approximations to
%  the integral form 0 to 1 of cos(alpha x^2) dx.

format short e

%  First compute integral for small alpha -- alpha = .1,.2,.3 --
%  and compare with Taylor series approximation. 

n = 100;
for alpha=.1:.1:.3,
  ftrap = trap(alpha,n);
  fsimp = simp(alpha,n);
  ftaylor = 1 - alpha^2/10 + alpha^4/216 - alpha^6/9360;
  alpha, ftrap, fsimp, ftaylor,
  input(' Hit return to continue ');
end;

%  Now do convergence study for alpha=1.

alpha = 1;

%  Use trapezoidal rule and Simpson's rule for n=2,4,8,...,128.

for i=1:7,
  n = 2^i;
  ftrap(i) = trap(alpha,n);
  fsimp(i) = simp(alpha,n);
end;

%%  Compute "exact" solution by using Simpson's rule with lots of points.
%%
%%n = 10*n;
%%fexact = simp(alpha,n);

%  Above method is ok, but let's do something more clever.
%  Compute "exact" solution by using Simpson's rule with Richardson
%  extrapolation.

fsimp(8) = simp(alpha,2*n);
fexact = (16*fsimp(8)-fsimp(7))/15;

%  Compute error in each approximation and ratios of successive errors.
%  Ratios should be about 4 for trapezoidal rule, about 16 for Simpson's
%  rule.

for i=1:7,
  errtrap(i) = ftrap(i)-fexact;
  errsimp(i) = fsimp(i)-fexact;
end;
for i=1:6,
  ratiotrap(i) = abs(errtrap(i)/errtrap(i+1));
  ratiosimp(i) = abs(errsimp(i)/errsimp(i+1));
end;

%  Print out errors and ratios.

errtrap, ratiotrap, 
errsimp, ratiosimp,