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.