Esercitazioni di MATLAB: Lezione del 7 Aprile 2008

Metodo delle secanti

Una possibile implementazione del metodo delle secanti:

File secanti.m

:open:

%
% Metodo delle secanti
%
% Cerco lo zero della funzione fun(x) a partire
% da due approssimazioni distinte dello zero x0 ed x1.
% epsi e' la tolleranza sotto la quale le iterate vengono
% interrotte.
%
% In uscita viene resttuito il vettore con tutte
% le iterate calcolate.
%
function X = secanti(fun,x0,x1,epsi)
  X  = [x0;x1] ;
  f0 = feval(fun,x0) ;
  f1 = feval(fun,x1) ;
  for i=1:100
    % formula di avanzamento
    df   = (f1-f0)/(X(end)-X(end-1)) ;
    xnew = X(end) - f1/df ;

    % aggiorna i vettori
    X  = [X ; xnew ] ;
    f0 = f1 ;
    f1 = feval(fun,xnew) ;
    
    % controlla se interrompere il ciclo
    if abs(f1) < epsi ; break ; end ;
  end
end

Uno script per plottare bene le iterate:

File piter.m

:open:

%
% Funzione per la stampa delle iterate per il metodo delle
% secanti
%
%  fun = funzione alla quale e' stato applicato
%        il metodo delle secanti
%  X   = vettore con la storia delle iterate
%
%  xe  = soluzione esatta
%
function piter(fun,X,xe)

  % calcolo ingombro iterate per coordinata X
  xmin = min(X) ;
  xmax = max(X) ;
  len  = xmax-xmin ;
  xmin = xmin - len/10 ;
  xmax = xmax + len/10 ;
  
  % vettori per il plot della funzione
  XX   = [xmin:len/1000:xmax] ;
  YY   = feval(fun,XX) ;

  % calcolo ingombro iterate per coordinata Y
  ymin = min(YY) ;
  ymax = max(YY) ;
  
  plot(XX,YY,'LineWidth',2);
  hold on;
  
  % assi centrati sullo zero
  plot([xe;xe],[ymin;ymax],'k-','LineWidth',1);
  plot([xmin;xmax],[0;0],'k-','LineWidth',1);
  
  for i=1:length(X)
    x = X(i) ;
    % linea verticale
    plot([x;x],[0;feval(fun,x)],'r--','LineWidth',2);
    if i>2
      x2 = X(i-2) ;
      % secante
      plot([x;x2],[0;feval(fun,x2)],'c-');      
    end
  end
  axis on ;
  hold off;
end

Uno script per testare tutto:

File test.m

:open:

%
% Funzione per testare il medoto delle secanti
%
%
function test
  
  % parametri per il problema
  x0   = 3 ;
  x1   = 1 ;
  epsi = 1E-15 ;
  
  % chiama il solutore
  X    = secanti(@fun,x0,x1,epsi) ;

  % soluzione esatta
  xe = 0;
  
  % differenza con la soluzione esatta
  ERR = X-xe

  % stima ordine di convergenza
  ORDER = log(ERR(3:end)./ERR(2:end-1)) ./ log(ERR(2:end-1)./ERR(1:end-2)) ;

  % stampa stima degli ordini
  disp(ORDER) ;
  
  % plottaggio delle iterate a video
  piter(@fun,X,xe) ;

end

function res = fun(x)
  res = x.*exp(x)./(1+x.^2) ;
  %res = (x-3).^3+(x-3)*0.1 ;
end

Tutti i file MATLAB in un unico archivio

lez3.zip

---