Analisi di stabilità lineare o assoluta stabilità

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

Considero il problema lineare:

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

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

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

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

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

ed è tale che

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

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

se $\mathop {Re} (\lambda )<0$$\mathop {Re} (\lambda )<0$.
Verifico la seguente definizione:

Definizione 7.7

Un metodo numerico è assolutamente stabile se

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

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

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

 

Definizione 7.8

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 \}$$

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

 

Suppongo che esistano due costanti $\mu _{min},\mu _{max}$$\mu _{min},\mu _{max}$ tali che $-\mu _{min}<f_{y}<-\mu _{max}$$-\mu _{min}<f_{y}<-\mu _{max}$, e pongo $\lambda =-\mu _{max}$$\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}$$

   $$y_{n+1}=y_{n}+hf_{n}$$
   se $f=\lambda y$$f=\lambda y$ segue che
   $$y_{n+1}=y_{n}+\lambda hy_{n}=(1+\lambda h)y_{n}=(1+\lambda h)^{n}y_{0}$$

   $$y_{n+1}=y_{n}+\lambda hy_{n}=(1+\lambda h)y_{n}=(1+\lambda h)^{n}y_{0}$$
   e per $n\to \infty$$n\to \infty$ la quantità tende a 0 se $|1+\lambda h|<1$$|1+\lambda h|<1$.Il metodo di eulero esplicito ha come condizione di assoluta stabilità $|1+\lambda h|<1$$|1+\lambda h|<1$.Nel piano di Argand-Gauss la regione di assoluta stabilità è un cerchio con raggio 1 e centro $1$$1$.$\lambda h$$\lambda h$ deve stare nel cerchio, altrimenti asintoticamente la soluzione numerica non si comporta come la soluzione esatta.In $\mathbb {R}$$\mathbb {R}$ la condizione diventa $h<-2/\lambda$$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}$$

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

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

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

   $$y_{n+1}={\frac {y_{0}}{(1-\lambda h)^{n}}}$$
   e impongo la condizione $|{\frac {1}{1-\lambda h}}|<1$$|{\frac {1}{1-\lambda h}}|<1$ che implica $|1-\lambda h|>1$$|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 considero il metodo di Frank-Nicholson,
   $$y_{n+1}=y_{n}+h/2(f_{n}+f_{n+1})$$

   $$y_{n+1}=y_{n}+h/2(f_{n}+f_{n+1})$$
   segue che, se $f=\lambda y$$f=\lambda y$, si ha
   $$y_{n+1}=y_{n}+h/2\lambda y_{n}+h/2\lambda y_{n+1}$$

   $$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}$$

   $$(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}$$

   $$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$$

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

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

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

   $$2h\lambda <0$$

Trasformazione di sistemi lineari

Considero il problema di Cauchy

$${\begin{cases}y'=F(t,y)\\y(t_{0})=y_{0}\end{cases}}$$

$${\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}$$F \colon \mathbb R \times \mathbb R^n \to \mathbb R^n$. E' un sistema di equazioni differenziali con

$$y'=Ay$$

$$y' = A y$$

Assumiamo che $A=S^{-1}\lambda S$$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)=\sum _{j=1}^{n}c_{j}e^{\lambda _{j}t}b_{j}$$

$y(t)\to 0\iff t\to \infty$$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'$$

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

$$y'=Sz'$$

$$y' = S z'$$

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

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

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

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

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

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

$$z_{i}'=\lambda _{i}z_{i}$$

$$z_{i}'=\lambda _{i}z_{i}$$

e considero il problema test su ciascuna delle equazioni. Scelgo come passo $h$$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}}$$

$${\begin{cases}y'=A\phi (x)\end{cases}}$$

con $\phi \in \mathbb {R} ^{n}$$\phi \in \mathbb {R} ^{n}$. Supponiamo che $A$$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)$$

$$y(t)=\sum _{j=1}^{n}c_{j}e^{\lambda _{j}t}v_{j}+\psi (t)$$

dove $\psi (t)$$\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}$$\lambda _{m},\lambda _{M}$ due autovalori della matrice $A$$A$ tali che valga la relazione:

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

$$|\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})|}}$$

$${\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$$\mathop {Re} (\lambda _{m})\sim 0$, anche se $\mathop {Re} (\lambda _{M})$$\mathop {Re} (\lambda _{M})$ non è molto grande, ho un problema stiff per definizione. I metodi utilizzati per risolverli sono i Runge-Kutta espliciti.