%
%  macm 202 -- 27 dec 02 -- djm
%
%  w1simp.m:  simpson's rule demo
%    - to run, type "w1simp" at the matlab prompt
%

%  clear workspace & graphics ("help clear", etc)
%  and initialize plot

clear;  close all;  disp(' ')
figure(1);  clf;  hold on

%  loop for 5 integrations
for m = 1:5;
    N = (2^m)*4;  L = sqrt(3);
    
%  interval discretization
    dt = L/N;  t = 0:dt:L;
    
%  simpson coefficent vector
    si = 1 - ((-1).^(0:1:N))/3;  si(1) = 1/3;  si(end) = 1/3;

%  simpson's rule
    integral = sum(w1simp1(t) .* si) * dt;
    
    error(m) = pi-integral;
    disp(['simp_' num2str(N) ': ' num2str(integral,8), ...
            '  ' num2str(error(m),4)])
    
%  plot accuracy
    plot(log10(4*(2.^m)),log10(error(m)),'x')
end
  
%  fine-tune error analysis plot

figure(1);  disp(' ')

xlabel('\bf log_{10} N    (N+1 = # of function evaluations)')
ylabel('\bf log_{10} |error| =  -(# digits accuracy)')
title('\bf log-log error plot')
plot([1.0 2.0],[-5.6 -9.6],'k')
text(1.2,-7.6,'\bf slope = -4')

axis([0.8 2.4 -13 -4])

%  numerical quadrature for decreasing tolerance (help quad)

[qintegral qN] = quad(@w1simp1,0,L,1e-6);
qerror = abs(pi-qintegral);
disp(['quad_' num2str(qN-1) ': ' num2str(qintegral,8), ...
            '  ' num2str(qerror,4)])
plot(log10(qN-1),log10(qerror),'ro')

[qintegral qN] = quad(@w1simp1,0,L,1e-8);
qerror = abs(pi-qintegral);
disp(['quad_' num2str(qN-1) ': ' num2str(qintegral,8), ...
            '  ' num2str(qerror,4)])
plot(log10(qN-1),log10(qerror),'ro')

[qintegral qN] = quad(@w1simp1,0,L,1e-10);
qerror = abs(pi-qintegral);
disp(['quad_' num2str(qN-1) ': ' num2str(qintegral,8), ...
            '  ' num2str(qerror,4)])
plot(log10(qN-1),log10(qerror),'ro')

% numerical quadrature using "quadl"

disp(' ')
[qintegral qN] = quadl(@w1simp1,0,L,1e-6);
qerror = abs(pi-qintegral);
disp(['quadl_' num2str(qN-1) ': ' num2str(qintegral,8), ...
            '  ' num2str(qerror,4)])
plot(log10(qN-1),log10(qerror),'k*')