Metodi numerici
per equazioni differenziali ordinarie - laboratorio
Docente:
Beatrice Meini
Lezione 6: Regioni di stabilita' assoluta; un esempio di metodo
predittore/correttore; confronto tra i vari
risolutori ode** di Matlab/Octave.
La
function stability_regionRK
disegna la frontiera delle regioni di stabilita' assoluta di vari metodi
espliciti di
Runge-Kutta. Osservare che le regioni sono piu' ampie quando il numero di stadi
cresce.
La function stability_regionAdams
disegna la frontiera delle regioni di stabilita' assoluta dei metodi
Adams-Moulton. Osservare che le regioni hanno intersezione sempre piu' piccola
con il semipiano complesso sinistro, mano a mano che il numero di passi k
aumenta.
Disegnare la frontiera delle regioni di stabilita' assoluta dei metodi
Adams-Bashforth (cercare ad esempio su Wikipedia i coefficienti).
- La function [y,t] = abm4(f,a,b,ya,n)
implementa un metodo predittore correttore di Adams-Bashforth-Moulton, di
ordine 4. Il metodo usa il metodo Adams-Bashforth a 4 passi come predittore,
il metodo di Adams-Moulton a 3 passi come correttore, e
un metodo di Runge-Kutta di ordine 4 come starter.
In input: f definisce l'equazione differenziale, a e b sono gli estremi
dell'intervallo d'integrazione, ya la condizione iniziale, e n il numero di
discretizzazioni dell'intervallo.
Applicare il metodo a problemi visti nelle lezioni precedenti. Come si comporta
il metodo applicato a problemi stiff?
- Il pendolo sferico:
Il moto di una particella di massa m, libera di muoversi all'interno di una sfera, e
soggetta solo alla forza di gravita', puo' essere descritto da una equazione
differenziale ordinaria di ordine 2. L'equazione e' ad esempio descritta nell'Example 2 di
questo documento.
La function y=fvinc(t,x), tratta dal libro
"Scientific Computing with MATLAB and Octave", Quarteroni, Saleri, Gervaso,
definisce l'equazione differenziale del primo ordine, che descrive il moto
della particella. Le prime 3 componenti di y rappresentano la derivata prima
della posizione in coordinate cartesiane, mentre le seconde 3 componenti sono
le rispettive derivate seconde; analogamente, le prime 3 componenti di x
rappresentano la posizione in coordinate cartesiane, mentre le seconde 3
componenti sono le rispettive derivate.
I vari risolutori ode** di Matlab/Octave si comportano in modo molto diverso
applicati a questo problema.
Risolvere il problema ai valori iniziali, con valore iniziale x0=[0, 1, 0, .8,
0, 1.2]; nell'intervallo [0 25], con vari ode** (ad esempio provare ode23 e ode45). Disegnare la traiettoria
della particella con il comando plot3; ad esempio, se x contiene l'output di
ode**, dare il comando plot3( x(:,1), x(:,2), x(:,3)). Cosa si osserva?
Provare a modificare le tolleranze.