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

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

Consideriamo il problema lineare:

${\begin{cases}y'=\lambda y\\y(0)=1\end{cases}}$

con $\lambda \in \mathbb {C}$. La soluzione esatta è

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

ed è tale che

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

se $\mathop {Re} (\lambda )<0$.

Verifichiamo la seguente definizione:

Definizione

Un metodo numerico è assolutamente stabile se

$|y_{n+1}|\to 0\iff n\to \infty$

(quindi se la soluzione numerica si comporta come la soluzione esatta).

Definizione

Si definisce la regione di assoluta stabilità

${\mathcal {A}}=\{z=\lambda h\in \mathbb {C} \,t.c.\,|y_{n+1}|\to 0\iff n\to \infty \}$

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}+hf_{n}$

se $f=\lambda y$ segue che

$y_{n+1}=y_{n}+\lambda hy_{n}=(1+\lambda h)y_{n}=(1+\lambda h)^{n}y_{0}$

e per $n\to \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 $\mathbb {R}$ la condizione diventa $h<-2/\lambda$.

2. considerando invece il metodo di Eulero implicito:

$y_{n+1}=y_{n}+hf_{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{cases}y'=F(t,y)\\y(t_{0})=y_{0}\end{cases}}$

con $F\colon \mathbb {R} \times \mathbb {R} ^{n}\to \mathbb {R} ^{n}$. È un sistema di equazioni differenziali con

$y'=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\iff t\to \infty$, e si richiede lo stesso comportamento per la soluzione numerica.

Allora si impone:

$z=S^{-1}y\,\longrightarrow \,z'=S^{-1}y'$

$y'=Sz'$

$y'=Ay\longrightarrow Sz'=ASz$

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

$z'=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}'=\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{cases}y'=A\phi (x)\end{cases}}$

con $\phi \in \mathbb {R} ^{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:

$|\mathop {Re} (\lambda _{M})|\geq |\mathop {Re} (\lambda _{j})|\geq \mathop {Re} (|\lambda _{m}|)$

Si definisce la stiffness ratio nel seguente modo:

${\frac {|\mathop {Re} (\lambda _{M})|}{|\mathop {Re} (\lambda _{m})|}}$

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

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