Esercitazioni di MATLAB: Lezione del 5 Maggio 2008

Controesempio di runge

Uno script per plottare l’esempio di Runge con nodi equispaziati:

File runge.m

:open:

%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.4 (Enrico Bertolazzi)
%
% Esempio di Runge per verificare la NON CONVERGENZA
% della intepolazione polinomiale.
%
function runge

  % metto i parametri per la interpolazione polinomiale in un cell array. 
  % I numeri sono il grado del polinomio interpolante, 
  % le stringhe sono la formattazione per il comanfo plot.
  N = { {4,'-k'}, ...
        {8,'-g'}, ...
        {12,'-m'} } ;
        % sopra il grado 12 ployfit da un avvertimento perche' non riesce
        % a calcolare "bene" il polinomio interpolante a causa del 
        % cattivo condizionamento del problema

  % x contiene le ascisse per il plottaggio delle curve e dei polinomi interpolanti
  x = -5:0.01:5;

  hold on ; % da adesso in poi il comando plot scrive sulla stessa finestra
  axis([-5 5 -1 1.1]) ; % dimensiono l'area di plottaggio
  plot( x, rfun(x), '-r', 'LineWidth', 2 ) ; % disegna la curva y = 1/(1+x^2)
  for k=1:length(N)
    % N{k}{1} = grado del polinomio interpolante
    % N{k}{2} = stringa formattazione per plot
    xp = -5:10/N{k}{1}:5 ; % costruisco vettore delle x con i punti di interpolazione equispaziati
    p  = polyfit( xp, rfun(xp), N{k}{1} ) ; % costruisco il polinomio interpolante,
                                            % fai doc polyfit per avere piu' informazioni
    plot( x, polyval(p, x), N{k}{2} ); % plottaggio del polinomio interpolante,
                                       % fai doc polyval per avere piu' informazioni
  end
  hold off ; % da adesso in poi il comando plot NON scrive sulla stessa finestra

end

function y = rfun( x )
  y = 1 ./ ( 1 + x .^ 2 ) ;
end

Uno script per plottare l’esempio di Runge con nodi di Chebishev:

File chebyshev.m

:open:

%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.4 (Enrico Bertolazzi)
%
% Esempio di Runge con nodi di Chebyshev.
%
function chebyshev

  % metto i parametri per la interpolazione polinomiale in un cell array. 
  % I numeri sono il grado del polinomio interpolante, 
  % le stringhe sono la formattazione per il comanfo plot.
  N = { {4,'-k'}, ...
        {8,'-g'}, ...
        {12,'-m'}, ...
        {16,'-b'} } ;
        % sopra il grado 12 ployfit da un avvertimento perche' non riesce
        % a calcolare "bene" il polinomio interpolante a causa del 
        % cattivo condizionamento del problema

  % x contiene le ascisse per il plottaggio delle curve e dei polinomi interpolanti
  x = -5:0.01:5;

  hold on ; % da adesso in poi il comando plot scrive sulla stessa finestra
  axis([-5 5 -0.2 1]) ; % dimensiono l'area di plottaggio
  plot( x, rfun(x), '-r', 'LineWidth', 2 ) ; % disegna la curva y = 1/(1+x^2)
  for k=1:length(N)
    % N{k}{1} = grado del polinomio interpolante
    % N{k}{2} = stringa formattazione per plot
    xp = 5*cheb( N{k}{1} ) % costruisco vettore delle x con i punti di interpolazione equispaziati
    p  = polyfit( xp, rfun(xp), N{k}{1} ) ; % costruisco il polinomio interpolante,
                                            % fai doc polyfit per avere piu' informazioni
    plot( x, polyval(p, x), N{k}{2} ); % plottaggio del polinomio interpolante,
                                       % fai doc polyval per avere piu' informazioni
  end
  hold off ; % da adesso in poi il comando plot NON scrive sulla stessa finestra

end

% fuzione di Runge
function y = rfun( x )
  y = 1 ./ ( 1 + x .^ 2 ) ;
end

% nodi di Chebyshev
function x = cheb( n )
  x = cos( [0:1/n:1] .* pi ) ;
end

Integrazione numerica

La funzione di Runge:

File rfun.m

:open:

%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.4 (Enrico Bertolazzi)
%
% Funzione di Runge
%
function y = rfun( x )
  y = 1 ./ ( 1 + x .^ 2 ) ;
end

Esercizio proposto dal prof. Ana Alonso

File alonso.m

:open:

%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.4 (Enrico Bertolazzi)
%
% Esercizio proposto dal prof. Ana Alonso
%
function alonso( N, a, b, fun )
  
  % nodi di interpolazione equispaziati
  xequ = [a:(b-a)/N:b] ;

  % nodi di interpolazione di Chebichev
  xcheb = (a+b)/2 + ((b-a)/2) * cheb(N) ;

  % polinomio interpolante con nodi equispaziati
  pequ  = polyfit( xequ, feval( fun, xequ), N ) ;

  % polinomio interpolante con nodi di Chebichev 
  pcheb  = polyfit( xcheb, feval( fun, xcheb), N ) ;

  % Primitiva del polinomio interpolante con nodi equispaziati
  Ipequ = polyint( pequ ) ;

  % Primitiva del polinomio interpolante con nodi di Chebichev
  Ipcheb = polyint( pcheb ) ;
  
  % calcolo integrale del polinomio pequ nell'intervallo [a,b]
  Iequ = polyval( Ipequ, b) - polyval( Ipequ, a) ;

  % calcolo integrale del polinomio pcheb nell'intervallo [a,b]
  Icheb = polyval( Ipcheb, b) - polyval( Ipcheb, a) ;

  % calcolo integrale della funzione fun con il comando quad
  Iquad = quad( fun, a, b ) ;
  
  % stampa dei risultati
  disp( sprintf( 'Integrale con nodi equispaziati    %f', Iequ) ) ;
  disp( sprintf( 'Integrale con nodi di Chebyshev    %f', Icheb) ) ;
  disp( sprintf( 'Integrale con funzione MATLAB quad %f', Iquad) ) ;
  disp( sprintf( 'Errore con nodi equispaziati       %f', Iequ-Iquad) ) ;
  disp( sprintf( 'Errore con nodi di Chebyshev       %f', Icheb-Iquad) ) ;

end

% nodi di Chebyshev nell'intervallo [-1,1]
function x = cheb( n )
  x = cos( [0:1/n:1] .* pi ) ;
end

Integrazione adattativa con il metodo dei trapezi:

File trapezi.m

:open:

%
% 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

Tutti i file MATLAB in un unico archivio

lez4.zip

---