% Programming and Numerical Methods
% MSc in Computational and Applied Physics

% Code 15: Comparing convergence of different quadrature formulae
% Midpoint, Trapezoidal, Cavalieri-Simpson , Newton-Cotes and Gauss-Legendre
% A. Meseguer (2007)

% Comparing convergence of different quadrature formulae
  clear all 
  close all 

% Selection interval [a,b]
  a = 0 ; b = 2*pi ;

  mp = [] ; t = [] ; cs = [] ; nc = [] ; gl = [] ;
  errmp = [] ; errt = [] ; errcs = [] ; errnc = [] ; errgl = [] ; 
  h = [] ;

  for m = [4:4:256]
   h = [h (b-a)/m] ;
   mp = [mp mpcf(a,b,m,@myfun)] ;
   t  = [t tcf(a,b,m,@myfun)] ;
   cs = [cs cscf(a,b,m,@myfun)] ;  
   nc = [nc cnc4f(a,b,m,@myfun)] ;
   gl = [gl glcf(a,b,m,@myfun)] ;
  end

% Following Quarteroni's example:
  exact = (3*(exp(-2*pi)-1) - 10*pi*exp(-2*pi))/25 ;
  errmp = abs(mp(1:end-1)-exact) ;
  errt = abs(t(1:end-1)-exact) ;
  errcs = abs(cs(1:end-1)-exact) ;
  errnc = abs(nc(1:end-1)-exact) ;
  errgl = abs(gl(1:end-1)-exact) ;

%  Otherwise, when we do NOT have the exact result
%  errmp = abs(mp(1:end-1)-mp(end)) ;
%  errt =  abs(t(1:end-1)-t(end)) ;
%  errcs = abs(cs(1:end-1)-cs(end)) ;
%  errnc = abs(nc(1:end-1)-nc(end)) ;
%  errgl = abs(gl(1:end-1)-gl(end)) ;

  h = h(1:end-1) ;
  loglog(h,errmp,'-or') ; hold on % Midpoint in red
  loglog(h,errt,'-ob') ;          % Trapezoidal in blue
  loglog(h,errcs,'-ok') ;         % Cavalieri-Simpson in black
  loglog(h,errnc,'-og') ;         % Newton-Cotes in green 
  loglog(h,errgl,'om') ;loglog(h,errgl,'-m') % Gauss-Legendre in magenta
  grid on