Lezione del 8 marzo 2007

Metodi di Runge-Kutta

Implementazione MATLAB di un generico metodo di Runge-Kutta esplicito:

File RK.m

:open:

%
% Metodo di Runge-Kutta 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]
%
%  Tableau struttura con campi
%   ns     = numero di stadi per il metodo di Runge Kutta
%   a(:,:) = matrice
%   b(:)   = vettore
%   c(:)   = vettore
%
% 
% Uscita
%  x(:)  = vettore soluzione delle 'x'
%  y(:)  = vettore soluzione delle 'y'
%
function [x,y]=RK(f,x0,y0,xf,N,TABLEAU)
  x(1) = x0 ;
  y(1) = y0 ;
  h    = (xf-x0)/N ;
  for k=1:N
    % griglia spazio
    x(k+1) = x(k) + h;
    for s=1:TABLEAU.ns
      % calcolo del coefficente Ks
      ys = y(k) ;
      for j=1:s-1
        ys = ys + TABLEAU.a(s,j)*K(j,:) ;
      end ;
      % xs = x(k) + c(s)*h
      % ys = y(k) + somma( a(s,j)* Kj)
      xs     = x(k)+TABLEAU.c(s)*h ;
      K(s,:) = h*feval(f,xs,ys) ;
    end;
    % y(k+1) = y(k) + somma( b(s) * Ks)
    y(k+1,:) = y(k,:) ;
    for s=1:TABLEAU.ns
      y(k+1,:) = y(k+1,:) + TABLEAU.b(s) * K(s,:) ;
    end ;
  end;
end

il codice ha bisogno del TABLEAU del metodo in ingresso.

Generazione del TABLEAU per un Runge-Kutta di ordine 3:

File RK3.m

:open:

function TABLEAU=RK3()
  TABLEAU.ns = 3 ;
  TABLEAU.c  = [0 1/3 2/3] ;
  TABLEAU.b  = [1/4 0 3/4] ;
  TABLEAU.a  = [ 0   0   0 ; ...
                 1/3 0   0 ; ...
                 0   2/3 0 ] ;
end

Generazione del TABLEAU per un Runge-Kutta di ordine 4:

File RK4.m

:open:

%
% Tableau di Butcher per Runge-Kutta
% al quarto ordine esplicito (classico)
%
function TABLEAU=RK4()
  TABLEAU.ns = 4 ;
  TABLEAU.c  = [0 1/2 1/2 1] ;
  TABLEAU.b  = [1/6 1/3 1/3 1/6] ;
  TABLEAU.a  = [ 0   0   0   0 ; ...
                 1/2 0   0   0 ; ...
                 0   1/2 0   0 ; ...
                 0   0   1   0 ] ;
end

Generazione del TABLEAU per il metodo di Heun:

File RKHeun.m

:open:

function TABLEAU=RKHeun()
  TABLEAU.ns = 2 ;
  TABLEAU.c  = [0,1] ;
  TABLEAU.b  = [1/2,1/2] ;
  TABLEAU.a  = [ 0 0 ; ...
                 1 0 ] ;
end

---