Calcolo numerico/Analisi di stabilità lineare o assoluta stabilità

Nella pratica, $h\to ̸0$, quindi la 0-stabilità non basta.

Consideriamo il problema lineare:

$\{\begin{array}{c}{y}^{\prime }=\lambda y\\ y(0)=1\end{array}$

con $\lambda \in ℂ$. La soluzione esatta è

$y(t)=e^{\lambda t}$

ed è tale che

$\underset{t\to \mathrm{\infty }}{lim}|y(t)|=0$

se $Re(\lambda )<0$.

Verifichiamo la seguente definizione:

Supponiamo che esistano due costanti ${\mu }_{min},{\mu }_{max}$ tali che $-{\mu }_{min}<f_y<-{\mu }_{max}$, e poniamo $\lambda =-{\mu }_{max}$.

La regione di assoluta stabilità varia a seconda del metodo considerato, infatti:

1. considerando il metodo di Eulero esplicito:

$y_{n+1}=y_n+h{f}_n$

se $f=\lambda y$ segue che

$y_{n+1}=y_n+\lambda h{y}_n=(1+\lambda h)y_n=(1+\lambda h{)}^n{y}_0$

e per $n\to \mathrm{\infty }$ la quantità tende a 0 se $|1+\lambda h|<1$. Il metodo di Eulero esplicito ha come condizione di assoluta stabilità $|1+\lambda h|<1$. Nel piano di Argand-Gauss la regione di assoluta stabilità è un cerchio con raggio 1 e centro $1$. $\lambda h$ deve stare nel cerchio, altrimenti asintoticamente la soluzione numerica non si comporta come la soluzione esatta. In $ℝ$ la condizione diventa $h<-2/\lambda$.

2. considerando invece il metodo di Eulero implicito:

$y_{n+1}=y_n+h{f}_{n+1}=y_n+h\lambda y_{n+1}$

$(1-\lambda h)y_{n+1}=y_n$

$y_{n+1}=\frac{y_n}{1-\lambda h}$

e ricorsivamente

$y_{n+1}=\frac{y_0}{(1-\lambda h{)}^n}$

e imponiamo la condizione $|\frac{1}{1-\lambda h}|<1$ che implica $|1-\lambda h|>1$. La regione di assoluta stabilità in questo caso è ciò che sta fuori dal cerchio di raggio 1 e centro 1. Quindi i metodi di Eulero implicito non hanno nessuna condizione sulla scelta del passo dovuti allo schema numerico.

3. Se consideriamo il metodo di Frank-Nicholson,

$y_{n+1}=y_n+h/2(f_n+f_{n+1})$

segue che, se $f=\lambda y$, si ha

$y_{n+1}=y_n+h/2\lambda y_n+h/2\lambda y_{n+1}$

$(1-h/2\lambda )y_{n+1}=y_n+h/2\lambda y_n=(1+h/2\lambda )y_n=(1+h/2\lambda {)}^n{y}_0$

quindi

$y_{n+1}=\frac{(1+h/2\lambda {)}^n}{1-h/2\lambda }{y}_0$

e si ottiene la condizione

$|\frac{1+h/2\lambda }{1-h/2\lambda }|=|\frac{2+h\lambda }{2-h\lambda }|<1$

$|2+h\lambda |<2-h\lambda$

$h\lambda -2<2+h\lambda <2-h\lambda$

$2h\lambda <0$

Trasformazione di sistemi lineari

Consideriamo il problema di Cauchy

$\{\begin{array}{c}{y}^{\prime }=F(t,y)\\ y(t_0)=y_0\end{array}$

con $F:ℝ\times {ℝ}^n\to {ℝ}^n$. È un sistema di equazioni differenziali con

$y^{\prime }=Ay$

Assumiamo che $A=S^{-1}\lambda S$ sia diagonalizzabile.

La soluzione esatta è

$y(t)=\sum _{j=1}^n{c}_j{e}^{{\lambda }_{j}t}{b}_j$

$y(t)\to 0⟺t\to \mathrm{\infty }$, e si richiede lo stesso comportamento per la soluzione numerica.

Allora si impone:

$z=S^{-1}y⟶z^{\prime }=S^{-1}{y}^{\prime }$

$y^{\prime }=S{z}^{\prime }$

$y^{\prime }=Ay⟶S{z}^{\prime }=ASz$

quindi, moltiplicando a sinistra per $S^{-1}$

$z^{\prime }=S^{-1}ASz=\lambda z$

si ha un sistema di equazioni differenziali lineari ma con $\lambda$ diagonale, in modo da avere equazioni differenziali disaccoppiate della forma:

$z_{i^{\prime }}={\lambda }_i{z}_i$

e consideriamo il problema test su ciascuna delle equazioni. Scegliamo come passo $h$, che deve valere per tutte le equazioni, l'autovalore di modulo maggiore.

Problemi stiff

Supponiamo di avere un sistema di equazioni differenziali del tipo

$\{\begin{array}{c}{y}^{\prime }=A\varphi (x)\end{array}$

con $\varphi \in {ℝ}^n$. Supponiamo che $A$ sia diagonalizzabile, e consideriamo il caso in cui le parti reali dei suoi autovalori siano negative.

Allora la soluzione esatta è nella forma

$y(t)=\sum _{j=1}^n{c}_j{e}^{{\lambda }_{j}t}{v}_j+\psi (t)$

dove $\psi (t)$ è La soluzione particolare della non omogenea.

Considerando un metodo numerico a una regione di assoluta stabilità limitata, allora ci sono necessariamente delle restrizioni sul passo di integrazione. Più piccola è la regione di assoluta stabilità, minore è il passo da scegliere, e la scelta del passo è dovuta alla soluzione transitoria. Il passo va scelto tanto più piccolo quanto più in fretta la soluzione tende a 0.

Chiamiamo ${\lambda }_m,{\lambda }_M$ due autovalori della matrice $A$ tali che valga la relazione:

$|Re({\lambda }_M)|\ge |Re({\lambda }_j)|\ge Re(|{\lambda }_m|)$

Si definisce la stiffness ratio nel seguente modo:

$\frac{|Re({\lambda }_M)|}{|Re({\lambda }_m)|}$

e il problema è stiff se questa quantità è molto maggiore di 1.

Se $Re({\lambda }_m)\sim 0$, anche se $Re({\lambda }_M)$ non è molto grande, abbiamo un problema stiff per definizione. I metodi utilizzati per risolverli sono i Runge-Kutta espliciti.