Esercitazioni di MATLAB: Lezione del 16 Maggio 2008

Problema proposto

EsercizioLab5.pdf

Soluzione proposta

Implementazione del metodo di Eulero:

File eulero.m

:open:

%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.5 (Enrico Bertolazzi)
%
% Solutore di Eulero
%
function [t,y] = eulero( fun, t0, y0, T, N )
  h = T/N ;       % calcola passo avanzamento
  t = t0:h:t0+T ; % calcola vettore tempi discrtizzati
  y(1) = y0 ;     % condizione iniziale
  for k=1:N
    % avanzamento con la formula di eulero
    y(k+1) = y(k) + h*feval(fun, t(k), y(k) ) ; 
  end ;
end

Funzione per testare il solutore:

Filefun.m

:open:

%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.5 (Enrico Bertolazzi)
%
% Funzione per ODE
%
function res = fun( t, y )
  res = 2*y./t + t.^2 .* exp(t) ;
end

Soluzione esatta per la funzione proposta:

File esatta.m

:open:

%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.5 (Enrico Bertolazzi)
%
% Soluzione esatta per ODE
%
function res = esatta( t )
  res = t.^2 .* ( exp(t) - exp(1) ) ;
end

Programma per testare il solutore:

File checkOrder.m

:open:

%
% Codici MATLAB per il corso di Calcolo Numerico 2007/2008
% Esercitazione n.5 (Enrico Bertolazzi)
%
% Controllo errore
clear all

Ns = { 5, 10, 20, 40, 80, 160, 320, 640 } ;

hold on

for k=1:length(Ns)
  [t,y] = eulero( @fun, 1, 0, 2, Ns{k} ) ;
  plot(t,y) ;
  err(k) = norm(y-esatta(t),inf) ;
end

t=1:0.001:3 ;
plot(t, esatta(t),'-r', 'LineWidth', 2 ) ;


disp('controllo a posteriori dell''errore' ) ;
for k=2:length(Ns)
  disp( sprintf('Ordine stimato %f', log(err(k-1)/err(k))/log(2) ) ) ;
end


hold off

Tutti i file MATLAB in un unico archivio

lez5.zip

---