Lezione del 1 marzo 2007

Esempio con una ODE elementare

Implementazione MATLAB del metodo di Eulero esplicito per la equazione di un circuito RL:

File circuito.m

:open:

%
% Soluzione esatta e approssimata con il metodo di Eulero
% della equazione differenziale
% L di(t)/dt + R i(t) = V
% i(0) = i0
%
% Ingresso:
%  i0 = corrente iniziale
%  V  = tensione ai capi del circuito
%  R  = resistenza
%  L  = induttanza
%  T  = tempo finale
%  N  = numero di punti di griglia da calcolare
% restituisce i vettori
% t(:) = vettore con gli istanti temporali dove la soluzione è calcolata
% i(:) = andamento della corrente calcolato con il metodo di Eulero
% e(:) = andamento della corrente con soluzione esatta
%
function [t,i,e]=circuito(i0,V,R,L,T,N)
  t(1) = 0  ;
  i(1) = i0 ;
  e(1) = i0 ;
  h    = T/N ;
  for k=1:N
    % griglia temporale
    t(k+1) = t(k) + h;
    % soluzione approssimata col metodo di Eulero
    i(k+1) = i(k) + (V-R*i(k))*h/L;
    % soluzione esatta
    e(k+1) = i0*exp(-t(k+1)*R/L) ...
             +(V/R)*(1-exp(-t(k+1)*R/L)) ;
  end;
end

Programma driver per verificare il codice:

File test_circuito.m

:open:

i0 = 1  ;
V  = 10 ;
R  = 1  ;
L  = 2  ;
T  = 10 ;
N  = 30 ;
[t,i,e] = circuito(i0,V,R,L,T,N);
hold off ;
plot(t,i,'.-b');
hold on ;
plot(t,e,'-r');
legend('approssimata','esatta');

Esempio di output:

Metodo di Eulero esplicito

Implementazione MATLAB del metodo di Eulero esplicito:

File Eulero_esplicito.m

:open:

%
% Metodo di Collatz o Eulero Migliorato per sistemi di equazioni differenziali
% 
% Ingresso
%  f      = funzione per il sistema y'=f(x,y)
%  x0, y0 = dato iniziale del problema y(x0) = y0
%  xf     = estremo superiore dell'intervallo di calcolo [x0,xf]
%  N      = numero di intervalli in cui suddividere  [x0,xf]
% 
% Uscita
%  x(:)  = vettore soluzione delle 'x'
%  y(:)  = vettore soluzione delle 'y'
%
function [x,y]=Eulero_esplicito(f,x0,y0,X,N)
 x(1)   = x0;
 y(1,:) = y0;
 h      = X/N ;
 for k=1:N
   % griglia spazio
   x(k+1) = x(k) + h;
   % soluzione approssimata col metodo di Eulero
   y(k+1,:) = y(k,:) + h*feval(f,x(k),y(k,:));
 end;
end

Metodo basato sullo sviluppo di Taylor al secondo ordine

ODE per l’esempio considerato:

File test_1_yp.m

:open:

%
% Esempio di metodo basato sullo sviluppo di Taylor
%
% funzione f(x,y)
%
function yp=test_1_yp(x,y)
  yp = -20.*y+20.*sin(x)+cos(x);
end

Derivata prima della ODE:

File test_1_ypp.m

:open:

%
% Esempio di metodo basato sullo sviluppo di Taylor
%
% funzione
%
%   d f(x,y)
%   -------
%    d x
%
function ypp=test_1_ypp(x,y)
  f   = - 20.*y+20.*sin(x)+cos(x);
  fx  =   20.*cos(x)-sin(x);
  fy  = - 20.*onex(size(x)) ;
  ypp = fx+fy.*f ;
end

Programma MATLAB che implementa il metodo basato sullo sviluppo di Taylor:

File test_1_taylor.m

:open:

%
% Esempio di metodo basato sullo sviluppo di Taylor
%
% calcolo della soluzione approssimata
%
function [x,y]=test_1_taylor(x0,y0,X,N)
  x(1) = x0;
  y(1) = y0;
  h    = X/N ;
  for k=1:N
    % griglia spazio
    x(k+1) = x(k) + h;
    % soluzione approssimata col metodo basato su Taylor
    yp  = test_1_yp(x(k),y(k)) ;
    ypp = test_1_ypp(x(k),y(k)) ;
    y(k+1) = y(k) + h*yp + (h*h/2)*ypp ;
  end;
end

Soluzione esatta del problema:

File test_1_esatta.m

:open:

%
% Esempio di metodo basato sullo sviluppo di Taylor
%
% soluzione esatta
%
function y=test_1_esatta(x)
  y=exp(-20.*x)+sin(x);
end

Programma driver per verificare il codice:

File test_taylor.m

:open:

x0 = 0 ;
y0 = 1 ;
X  = 2 ;
N  = 30 ;
[x,y]=test_1_taylor(x0,y0,X,N) ;
hold off ;
plot(x,y,'.-b');
hold on ;
x = [x0:(X-x0)/200:X];
plot(x,test_1_esatta(x),'-r');
legend('approssimata','esatta');

Esempio di output:

Metodo al secondo ordine

Implementazione MATLAB del metodo di Collatz:

File Collatz.m

:open:

%
% Metodo di Collatz o Eulero Migliorato per sistemi di equazioni differenziali
% 
% Ingresso
%  f      = funzione per il sistema y'=f(x,y)
%  x0, y0 = dato iniziale del problema y(x0) = y0
%  xf     = estremo superiore dell'intervallo di calcolo [x0,xf]
%  N      = numero di intervalli in cui suddividere  [x0,xf]
% 
% Uscita
%  x(:)  = vettore soluzione delle 'x'
%  y(:)  = vettore soluzione delle 'y'
%
function [x,y]=Collatz(f,x0,y0,xf,N)
  x(1)   = x0;
  y(1,:) = y0;
  h      = (xf-x0)/N ;
  for k=1:N
    % griglia spazio
    x(k+1) = x(k) + h;
    % mezzo passo di Eulero
    % y(k+1/2) = y(k) + (h/2)*f(x(k),y(k))
	fk  = feval(f,x(k),y(k,:));
    ykh = y(k,:) + (h/2).*fk;
    % correzione di Collatz
    % y(k+1) = y(k) + h*f(x(k+1/2),y(k+1/2))
	fkh = feval(f,x(k) + h/2,ykh);
    y(k+1,:) = y(k,:) + h.*fkh ;
  end;
end

Implementazione MATLAB del metodo di Heun:

File Heun.m

:open:

%
% Metodo di Heun per sistemi di equazioni differenziali
% 
% Ingresso
%  f      = funzione per il sistema y'=f(x,y)
%  x0, y0 = dato iniziale del problema y(x0) = y0
%  xf     = estremo superiore dell'intervallo di calcolo [x0,xf]
%  N      = numero di intervalli in cui suddividere  [x0,xf]
% 
% Uscita
%  x(:)  = vettore soluzione delle 'x'
%  y(:)  = vettore soluzione delle 'y'
%
%
function [x,y]=Heun(f,x0,y0,xf,N)
  x(1)   = x0 ;
  y(1,:) = y0 ;
  h      = (xf-x0)/N ;
  for k=1:N
    % griglia spazio
    x(k+1) = x(k) + h;
    % passo di Eulero
    % ytilde(k+1) = y(k) + h*f(x(k),y(k))
	fk     = feval(f,x(k),y(k,:));
    ytilde = y(k,:) + h*fk;
    % correzione di Heun
    % y(k+1) = y(k) + (h/2)*(f(x(k),y(k))+f(x(k+1),ytilde(k+1)))
	fk1      = feval(f,x(k+1),ytilde);
    y(k+1,:) = y(k,:) + h*(fk+fk1)/2 ;
  end;
end

Programma driver per verificare i metodi di Collatz, Heun e Eulero:

File test_multiplo.m

:open:

x0 = 0 ;
y0 = 1 ;
X  = 2 ;
N  = 30 ;
[x,ye] = Eulero_esplicito (@test_1_yp,x0,y0,X,N) ;
[x,yc] = Collatz          (@test_1_yp,x0,y0,X,N) ;
[x,yh] = Heun             (@test_1_yp,x0,y0,X,N) ;
hold off ;
plot(x,ye,'.-b');
hold on ;
plot(x,yc,'.-k');
plot(x,yh,'-gs');
x = [x0:(X-x0)/200:X];
plot(x,test_1_esatta(x),'-r','LineWidth',2);
legend('Eulero',  ...
       'Collatz', ...
       'Heun',    ...
       'Esatta');

Esempio di output:

Esempio per un sistema di ODE

Definizione della ODE:

File cerchio_yp.m

:open:

%
% Equazioni per il sistema differenziale
%
% x'  = y
% y'  = -x
%
function Zp=cerchio_yp(t,Z)
  Zp(1) = Z(2);
  Zp(2) = -Z(1);
end

Programma driver per verificare i metodi di Collatz, Heun e Eulero sulla ODE:

File test_cerchio.m

:open:

x0 = 0  ;
y0 = 1  ;
X  = 10 ;
N  = 30 ;
[t,e] = Eulero_esplicito (@cerchio_yp,0,[x0,y0],X,N*2) ;
[t,c] = Collatz          (@cerchio_yp,0,[x0,y0],X*2,N) ;
[t,h] = Heun             (@cerchio_yp,0,[x0,y0],X*2,N) ;
hold off ;
plot(e(:,1),e(:,2),'.-b');
hold on ;
plot(c(:,1),c(:,2),'.-k');
plot(h(:,1),h(:,2),'-gs');
t = [0:0.01:2*pi];
x = cos(t) ;
y = sin(t) ;
plot(x,y,'-r','LineWidth',2);
legend('Eulero',  ...
       'Collatz', ...
       'Heun',    ...
       'Esatta');

Esempio di output:

---