Metodi numerici
per equazioni differenziali ordinarie - laboratorio
Docente:
Beatrice Meini
Lezione 4: Problemi stiff
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 metodi di Runge-Kutta impliciti, generalmente adatti a risolvere problemi stiff. In particolare Matlab ha il risolutore ode23s, che implementa un metodo di R-K con convergenza 2/3. Octave ha i risolutori odebwe, che implementa il metodo di Eulero implicito, ode2r e ode5r, che implementano metodi di Radau con diverso numero di stadi.
La function [y] = es_stiff1(t,x,lambda) e' la funzione che definisce l'equazione differenziale vista a lezione, che verifica l'intensita' di corrente in un circuito elettrico. L'IVP con condizione iniziale x(0)=0 e intervallo [0 5] e' stiff se lambda e' reale, negativo, e molto grande in modulo. Pet t "non vicino" a zero, la soluzione e' ben approssimata da cos(t). Confrontare ode45 con odebwe, ode2r, ode5r in Octave, oppure con ode23s in Matlab. In particolare confrontare le statistiche e disegnare le approssimazioni ottenute. Osservare, facendo degli zoom, che ode45 produce una approssimazione che oscilla e che richiede moltissimi punti di discretizzazione.
La function [y] = es_stiff2(t,x) e' la funzione che definisce una equazione differenziale. Fissiamo un numero reale positivo delta. L'IVP con condizione iniziale x(0) = delta e intervallo [ 0 2/delta] e' stiff se delta e' piccolo. L'intervallo "critico" in questo caso e' in un intorno del punto di mezzo dell'intervallo di integrazione. Anche in questo esempio confrontare i vari risolutori, seguendo le indicazioni all'interno del file.
La function [y] = es_stiff3(t,x) modellizza una reazione chimica tra 2 sostanze. L'IVP con condizione iniziale [1; 0] e intervallo [ 0 50] e' stiff. In questo esempio ode45 e ode78 di Octave, con i parametri di default, falliscono. Anche in questo esempio confrontare i vari risolutori, seguendo le indicazioni all'interno del file. In particolare utilizzare lo Jacobiano.
Le function rcd (dal libro MATLAB Guide, di D.J. Higham e N.J. Higham) e rcd_octave (la function rcd adattata a Octave) risolvono un IVP N-dimensionale stiff, passando tra le opzioni lo Jacobiano e la sua struttura di sparsita'. L'IVP e' descritto qui. Provare a modificare il risolutore, le opzioni e le dimensioni.