%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.4 (Enrico Bertolazzi)
%
% Funzione che calcola con preciszione stimata epsi l'integrale della funzione fun
%
% fun  = puntatore alla funzione da integrare
% a    = estremo sinistro di integrazione
% b    = estremo destro di integrazione
% epsi = tolleranza
%
function res = trapezi( fun, a, b, epsi ) 

  % alcuni controlli sei dati
  if b    < a     ; error( 'intervallo di integrazione errato b < a!' ) ; end ;
  if epsi < 1E-10 ; error( 'tolleranza troppo stretta (<1E-10) eo negativa' ) ; end ;

  % calcolo integrale con un solo trapezio
  n   = 1 ;
  res = trap( fun, a, b, n ) ;
  
  for k=1:10
    % salva l'integrale precedentemente calcolato
    res_old = res ;
    
    % raddoppia il numero di intervalli
    n       = n * 2 ;
    res     = trap( fun, a, b, n ) ; 
    
    % stima dell'errore con la formula di estrapolazione
    err     = (res - res_old) / 3 ;
    
    % se l'errore stimato e' sotto la tolleranza
    % interrompo le iterate
    if abs( err ) < epsi ; break ; end ;
  end
end

%
% Calcola l'integrale approssimato con la formula dei trapezi
%
function res = trap( fun, a, b, n )
  % nodi quadratura
  h     = (b-a)/n ;
  xnode = [a+h:h:b-h] ;
  % contributo nodi interni
  res   = h * sum(feval(fun, xnode)) ;
  % contributo estremi
  res   = res + (h/2)*( feval( fun, a ) + feval( fun, b ) ) ;
end