Metodi numerici per equazioni differenziali ordinarie - laboratorio
Docente: Beatrice Meini

Lezione 3: I metodi di Runge-Kutta implementati da Matlab/Octave

Attenzione: La versione 3.2 di Octave ha un package per risovere ODE; per utilizzare questa versione occorre dare il comando octave3.2.

• Matlab e Octave hanno dei risolutori di IVP che implementano vari metodi di Runge-Kutta. Matlab ha ode23 e ode45 che implementano 2 metodi di R-K embedded, con convergenza di ordine 2/3 e 4/5, rispettivamente. Octave, oltre a ode23 e ode45, ha ode54 e ode78; alle pagine 6-7 del manuale trovate una descrizione dei metodi implementati; in particolare ode78 ha ordine di convergenza 7/8, e ode54 (che ha ordine di convergenza 5/4) fornisce una migliore stima dell'errore locale di discretizzazione.
  
  L'uso elementare e' [t, x] = odeXY(@fname, tspan, x0), dove fname e' il nome della function che definisce l'equazione differenziale, tspan e' l'intervallo d'integrazione, e x0 e' il valore iniziale.
  
  Usi piu' avanzati si ottengono modificando le opzioni, mediante il comando odeset (si vedano le possibili opzioni con l'help di Matlab, o alla pagina 10 del manuale di Octave). Ad esempio il comando opzioni=odeset('AbsTol',1.e-8,'RelTol',1.e-7,'Stats','on'); modifica le tolleranze per il controllo dell'errore, e attiva la visualizzazione delle statistiche. Le opzioni vengono caricate dando il comando [t, x] = odeXY(@fname, tspan, x0, opzioni).
  
  In Octave il comando S=odeXY(@fname, tspan, x0, opzioni) genera un struct array, tale che la variabile S.x contiene la discretizzazione di tspan, variabile S.y contiene le approssimazioni della soluzione, e la variabile S.stats contiene le statistiche. In Matlab le statistiche vengono mostrate nella command window.

• Provare i vari risolutori odeXY che implementano metodi di RK, applicati agli IVP definiti dalla function [y] = f3(t,x) e dalla function [y] = moto(t,x) della Lezione 2. In particolare, osservare come viene discretizzato l'intervallo, modificando le tolleranze per il controllo dell'errore locale; questo di osserva facilmente aggiungendo l'opzione '+', oppure 'o', nel comando plot che traccia il grafico della soluzione. Osservare le statistiche, in particolare i fallimenti e il numero di valutazioni della funzione.

• La function [y] = f4(t,x,lambda) e' la funzione che definisce una equazione differenziale in cui i metodi odeXY falliscono, in quando l'errore globale e' molto piu' grande dell'errore locale. La soluzione dell'IVP con x0=0 e' la funzione sin(t). Ad esempio provare con lambda=50:
  
  opzioni=odeset('AbsTol',1.e-6,'RelTol',1.e-4,'Stats','on');
  S=ode45(@f4,[0 10], 0,opzioni,50);
  plot(S.x, S.y, S.x, sin(S.x))
  Cosa si osserva? Osservare le statistiche: i fallimenti sono trascurabili, dunque il risolutore pensa di aver fornito una buona approssimazione.
  Provare altri odeXY.

• Le opzioni permettono di interrompere l'integrazione dell'IVP quando si determina un particolare "evento", definito da una function. Un esempio e' il modello di una volpe che insegue un coniglio, seguendo una traiettoria definita da una equazione differenziale. Se la volpe raggiunge il coniglio voglio interrompere il calcolo della traiettoria, anche se l'intervallo di integrazione non e' stato completato. Infatti, quando la volpe raggiunge il coniglio l'equazione differenziale non e' piu' definita. Dettagli sul modello e sull'equazione differenziale si trovano qui. Il file volpe.tar contiene delle function (testate in Matlab, ma NON in octave, per cui la sintassi va modificata) per risolvere l'IVP con e senza il controllo dell'evento. In particolare, la function [tfox,yfox]=fox_rabbit1(k) non usa il controllo dell'evento e, se la volpe raggiunge il coniglio, il programma termina con un errore. La function [tfox,yfox]=fox_rabbit2(k) invece usa il controllo dell'evento e, se la volpe raggiunge il coniglio, il programma termina l'integrazione quando si verifica l'evento, senza errore.