31  Stabilità assoluta

Autore/Autrice

Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

L’accuratezza da sola non basta a decidere se un metodo numerico per ODE sia utilizzabile. Per problemi che contengono modi che decadono rapidamente un metodo può essere di ordine elevato e ciò nonostante richiedere un passo h piccolissimo solo per non esplodere. Questo capitolo sviluppa il concetto di stabilità assoluta e confronta le regioni di stabilità dei metodi visti nei capitoli precedenti.

ImportanteTre nozioni da non confondere

• Buona posizione (27.7): proprietà del problema continuo.
• Convergenza (Lemma 28.1): comportamento dell’errore per h\to 0 su un intervallo finito.
• Stabilità assoluta: comportamento della soluzione numerica per n\to\infty a passo h fissato.

Sono indipendenti: un metodo convergente può essere inutilizzabile su un problema stiff, e un metodo stabile può essere del tutto impreciso.

31.1 L’equazione test di Dahlquist

Il modello di riferimento è l’equazione test scalare

y'(t) = \lambda\, y(t), \qquad \lambda\in\mathbb{C} \tag{31.1}

la cui soluzione esatta avanza di un passo secondo

y(t_{n+1}) = e^{q}\, y(t_n), \qquad q = h\lambda

Se \mathrm{Re}(\lambda)<0 la soluzione esatta decade. È ragionevole chiedere che un buon metodo numerico riproduca questo decadimento per il valore q=h\lambda scelto.

Perché proprio questa equazione? Perché, come vedremo nella Sezione 31.5, la linearizzazione di un sistema attorno ad una traiettoria e la diagonalizzazione della matrice jacobiana riconducono il caso vettoriale ad una collezione di equazioni scalari della forma 31.1, una per ogni autovalore.

31.2 Funzione e regione di stabilità

Applicando un metodo a un passo alla 31.1 si ottiene sempre una relazione della forma

y_{n+1} = R(q)\, y_n \qquad\Longrightarrow\qquad y_n = R(q)^{n}\, y_0

dove R(q) è detta funzione di stabilità del metodo. La soluzione numerica resta limitata per n\to\infty se e solo se |R(q)|\leq 1.

Definizione 31.1 (Regione di stabilità assoluta) Dato un metodo a un passo con funzione di stabilità R(q), la sua regione di stabilità assoluta è

\mathcal{S} = \big\lbrace q\in\mathbb{C} ~:~ |R(q)| \leq 1 \big\rbrace \tag{31.2}

La frontiera di \mathcal{S} è descritta da |R(q)|=1. Per ogni metodo consistente vale

R(q) = 1 + q + \mathcal{O}(q^2)

cioè R riproduce localmente l’esponenziale e^q vicino all’origine; più precisamente, un metodo di ordine p ha

R(q) = \sum_{k=0}^{p} \dfrac{q^k}{k!} + \mathcal{O}(q^{p+1})

ConsiglioStabilità dal tableau di Butcher

Per un metodo di Runge-Kutta¹ con tableau (\bm{A},\bm{b},\bm{c}) a s stadi (Capitolo 30) applicato alla 31.1 si ha

\bm{Z} = y_n\bm{1} + q\bm{A}\bm{Z} \qquad\Longrightarrow\qquad \bm{Z} = y_n(\bm{I}-q\bm{A})^{-1}\bm{1}

e quindi

\boxed{\; R(q) = 1 + q\, \bm{b}^{T} (\bm{I} - q\bm{A})^{-1} \bm{1} \;} \tag{31.3}

Se il metodo è esplicito, \bm{A} è strettamente triangolare inferiore, quindi \bm{A}^{s}=\bm{0} e la serie di Neumann² (\bm{I}-q\bm{A})^{-1}=\sum_{k\geq 0}q^k\bm{A}^k è finita: R(q) è un polinomio di grado \leq s. Se il metodo è implicito, R(q) è una funzione razionale.

Questa è l’osservazione decisiva: un polinomio non costante è illimitato, quindi |R(q)|\leq 1 non può valere su un insieme illimitato. Ne segue che

nessun metodo di Runge-Kutta esplicito può avere regione di stabilità illimitata.

31.3 I metodi fondamentali

Eulero esplicito

y_{n+1} = y_n + h\lambda y_n = (1+q) y_n \qquad\Longrightarrow\qquad \boxed{\,R_{\mathrm{EE}}(q) = 1+q\,}

La condizione |1+q|\leq 1 descrive il disco chiuso di raggio 1 centrato in -1. Sull’asse reale negativo

-2 \leq q \leq 0

quindi per \lambda reale negativo il metodo richiede

h \leq \dfrac{2}{|\lambda|} \tag{31.4}

Si tratta di un vincolo di pura stabilità: può essere enormemente più severo del passo necessario a risolvere accuratamente la parte lenta della soluzione.

Runge-Kutta espliciti del secondo ordine

Tutti i metodi RK espliciti a due stadi di ordine due (punto medio, Heun³, Ralston⁴) hanno la stessa funzione di stabilità

\boxed{\,R_{\mathrm{RK2}}(q) = 1 + q + \dfrac{q^2}{2}\,}

I coefficienti dei tableau sono diversi, ma sull’equazione test i metodi coincidono. L’intervallo di stabilità sull’asse reale negativo è ancora [-2,0]: aumentare l’ordine da 1 a 2 non allarga il passo ammissibile per un modo reale negativo.

Runge-Kutta classico del quarto ordine

\boxed{\, R_{\mathrm{RK4}}(q) = 1 + q + \dfrac{q^2}{2} + \dfrac{q^3}{6} + \dfrac{q^4}{24} \,}

cioè il polinomio di Taylor⁵ di grado 4 di e^q. L’intervallo di stabilità reale è circa

-2.7853 \lesssim q \leq 0

un po’ più largo di Eulero⁶, ma pur sempre limitato.

Eulero implicito

y_{n+1} = y_n + q y_{n+1} \qquad\Longrightarrow\qquad \boxed{\,R_{\mathrm{EI}}(q) = \dfrac{1}{1-q}\,}

La condizione |R|\leq 1 diventa |1-q|\geq 1, cioè l’esterno del disco di raggio 1 centrato in +1. Tale insieme contiene tutto il semipiano sinistro chiuso. Inoltre

R_{\mathrm{EI}}(q) \longrightarrow 0 \qquad \text{per } |q|\to\infty,\ \mathrm{Re}(q)<0

quindi i modi molto rigidi vengono fortemente smorzati.

Metodo dei trapezi (Crank-Nicolson)

y_{n+1} = y_n + \dfrac{q}{2}\big( y_n + y_{n+1} \big) \qquad\Longrightarrow\qquad \boxed{\, R_{\mathrm{CN}}(q) = \dfrac{1 + q/2}{1 - q/2} \,}

Si verifica subito che |R_{\mathrm{CN}}(q)|\leq 1 se e solo se \mathrm{Re}(q)\leq 0: la regione di stabilità è esattamente il semipiano sinistro. Però

R_{\mathrm{CN}}(q) \longrightarrow -1 \qquad \text{per } q\to -\infty

I modi rigidi non vengono smorzati: sopravvivono numericamente come oscillazioni di segno alterno. È questa la differenza pratica fondamentale fra trapezi ed Eulero implicito.

31.4 A-stabilità e L-stabilità

Definizione 31.2 (Metodo A-stabile) Un metodo a un passo si dice A-stabile se

\big\lbrace q\in\mathbb{C} ~:~ \mathrm{Re}(q)\leq 0 \big\rbrace \subseteq \mathcal{S}

cioè se ogni modo stabile dell’equazione test resta numericamente stabile per ogni passo h>0.

Definizione 31.3 (Metodo L-stabile) Un metodo A-stabile si dice L-stabile se in più

R(q) \longrightarrow 0 \qquad \text{per } |q|\to\infty \text{ nel semipiano sinistro}

L’A-stabilità impedisce la crescita dei modi rigidi stabili; la L-stabilità li sopprime. Dalla discussione precedente:

• Eulero implicito è A-stabile e L-stabile;
• il metodo dei trapezi è A-stabile ma non L-stabile;
• nessun metodo di Runge-Kutta esplicito è A-stabile (31.3).

ImportanteBarriera di Dahlquist

Si può dimostrare (seconda barriera di Dahlquist⁷) che un metodo lineare multipasso A-stabile ha ordine al più 2, e che fra questi il metodo dei trapezi è quello con la costante d’errore minima. Per i metodi di Runge-Kutta impliciti invece non c’è barriera: i metodi di Gauss-Legendre⁸ (Capitolo 30) sono A-stabili a ogni ordine, i metodi di Radau⁹ IIA sono L-stabili. La dimostrazione si trova in [3:V].

31.5 Dal caso scalare ai sistemi

Si consideri il sistema lineare

\bm{y}' = \bm{A}\bm{y}

Per un metodo di Runge-Kutta la funzione di stabilità si estende formalmente all’argomento matriciale:

\bm{y}_{n+1} = R(h\bm{A})\, \bm{y}_n

Se \bm{A} è diagonalizzabile, \bm{A}=\bm{V}\bm{\Lambda}\bm{V}^{-1}, allora

R(h\bm{A}) = \bm{V}\, R(h\bm{\Lambda})\, \bm{V}^{-1}

e poiché R(h\bm{\Lambda}) è diagonale, il problema si disaccoppia in tante equazioni test scalari quanti sono gli autovalori. La condizione pratica è quindi

\boxed{\; h\lambda_j \in \mathcal{S} \quad \text{per ogni autovalore } \lambda_j \text{ di } \bm{A} \;} \tag{31.5}

AvvisoLimiti del criterio spettrale

La 31.5 è esatta per matrici normali. Per matrici fortemente non normali gli autovalori da soli possono non bastare: la crescita transitoria è governata dal numero di condizionamento di \bm{V} e non dallo spettro. I diagrammi di stabilità vanno quindi letti come diagnostica, non come teoria completa, soprattutto per problemi non lineari.

Problemi stiff

Si consideri

\bm{y}' = \bm{A}\bm{y}, \qquad \bm{A} = \begin{bmatrix} -1 & 0 \\ 0 & -1000 \end{bmatrix}

I due modi decadono come e^{-t} ed e^{-1000t}. Dopo un tempo brevissimo il secondo modo è numericamente nullo e la soluzione è governata interamente dal primo: per accuratezza basterebbe h\approx 10^{-1}. Ma Eulero esplicito, per la 31.4, richiede

h \leq \dfrac{2}{1000} = 2\cdot 10^{-3}

per non esplodere, cioè un passo 500 volte più piccolo di quello richiesto dall’accuratezza, e lo richiede per tutta la durata dell’integrazione. Un problema con questa caratteristica si dice stiff; una misura grossolana è il rapporto

\dfrac{\max_j |\mathrm{Re}(\lambda_j)|}{\min_j |\mathrm{Re}(\lambda_j)|}

detto rapporto di stiffness (qui vale 1000). Per questi problemi i metodi impliciti A-stabili non sono un lusso ma una necessità.

31.6 Diagrammi di stabilità

Le regioni che seguono sono calcolate direttamente dalla definizione 31.2, tracciando la curva di livello |R(q)|=1.

Figura 31.1: Regioni di stabilità assoluta di Eulero esplicito, RK2 e RK4 classico. La regione stabile è l’interno della curva.

Le regioni sono tutte limitate: aumentare l’ordine ne cambia forma ed estensione, ma non produce A-stabilità.

Figura 31.2: Frontiere di stabilità di Eulero implicito e del metodo dei trapezi. Entrambi contengono il semipiano sinistro; solo Eulero implicito è L-stabile.

AvvisoCome si legge il disegno

Per Eulero implicito la frontiera è la circonferenza |1-q|=1, ma la parte stabile è l’esterno del cerchio, non l’interno. Per il metodo dei trapezi la frontiera è l’asse immaginario e la parte stabile è tutto il semipiano sinistro. Una curva di livello, da sola, non dice da che parte si sta: va sempre letta insieme alla definizione 31.2.

Il grafico seguente mostra il confronto fra |R(q)| ed |e^q|=e^{q} sull’asse reale negativo: si vede dove ogni metodo perde il contatto con la soluzione esatta.

Figura 31.3: Modulo della funzione di stabilità sull’asse reale negativo confrontato con la soluzione esatta e^q. La soglia |R|=1 delimita la stabilità.

31.7 Stabilità e accuratezza

Il fatto che q appartenga a \mathcal{S} garantisce solo che la propagazione numerica resti limitata: non dice che l’approssimazione sia accurata. Eulero implicito è stabile per q negativo arbitrariamente grande, ma

R_{\mathrm{EI}}(q) = \dfrac{1}{1-q}

approssima male e^q a passo grossolano. La scelta pratica di h deve quindi soddisfare due condizioni separate:

1. stabilità: gli autovalori scalati h\lambda_j devono stare in \mathcal{S};
2. accuratezza: l’errore globale deve rispettare la tolleranza richiesta.

Per problemi non stiff è di norma l’accuratezza a determinare h; per problemi stiff un metodo esplicito è invece limitato dalla stabilità.

Problemi oscillanti e asse immaginario

Per problemi quasi conservativi gli autovalori stanno vicino all’asse immaginario e conta il comportamento di R(i\omega h). Si distinguono due errori:

• errore di dissipazione: |R(i\omega h)| \neq 1;
• errore di dispersione (o di fase): \arg R(i\omega h) \neq \omega h.

Un metodo può quindi essere assolutamente stabile e riprodurre comunque male le oscillazioni. Si noti che Eulero esplicito ha |R(i\theta)|=\sqrt{1+\theta^2}>1: nessun tratto dell’asse immaginario, tranne l’origine, appartiene alla sua regione di stabilità, per cui esso è inutilizzabile sui problemi puramente oscillanti.

Procedura pratica

1. Linearizzare localmente la ODE: \bm{y}' \approx \bm{J}_{\bm{f}}(t,\bm{y})\,\bm{y} + \cdots
2. Stimare gli autovalori \lambda_j dello jacobiano \bm{J}_{\bm{f}}.
3. Scalarli con il passo candidato: q_j = h\lambda_j.
4. Verificare che i q_j stiano in \mathcal{S}.
5. Imporre separatamente la tolleranza di accuratezza.

31.8 Quadro riassuntivo

Tabella 31.1: Confronto fra i metodi a un passo studiati.

Metodo	Ordine	Tipo	R(q)	Stabilità su \mathbb{R}^{-}	A-stab.	L-stab.
Eulero esplicito	1	esplicito	1+q	[-2,0]	no	no
RK2 (Heun, punto medio, Ralston)	2	esplicito	1+q+\dfrac{q^2}{2}	[-2,0]	no	no
RK4 classico	4	esplicito	\sum_{k=0}^{4}\dfrac{q^k}{k!}	\approx[-2.7853,\,0]	no	no
Eulero implicito	1	implicito	\dfrac{1}{1-q}	(-\infty,0]	sì	sì
Trapezi (Crank-Nicolson)	2	implicito	\dfrac{1+q/2}{1-q/2}	(-\infty,0]	sì	no
Punto medio implicito	2	implicito	\dfrac{1+q/2}{1-q/2}	(-\infty,0]	sì	no
Gauss-Legendre (s stadi)	2s	implicito	Padé (s,s)	(-\infty,0]	sì	no
Radau IIA (s stadi)	2s-1	implicito	Padé (s-1,s)	(-\infty,0]	sì	sì

I punti essenziali del capitolo:

• le regioni di stabilità dei metodi di Runge-Kutta espliciti sono sempre limitate, perché R è un polinomio;
• un metodo A-stabile contiene tutto il semipiano sinistro, e deve necessariamente essere implicito;
• la L-stabilità aggiunge lo smorzamento dei modi molto rigidi;
• il metodo dei trapezi è A-stabile ma non L-stabile;
• su un problema stiff il passo di un metodo esplicito è dettato dalla stabilità, non dall’accuratezza;
• stabilità e accuratezza sono requisiti indipendenti.

[1]

J. C. Butcher. 1987. The Numerical analysis of ordinary differential equations : Runge-Kutta and general linear methods (First ed.). Wiley, Chichester.

[2]

Frank Cameron. 2006. A Matlab package for automatically generating Runge-Kutta trees, order conditions, and truncation error coefficients. ACM Trans. Math. Softw. 32, 2 (giugno 2006), 274–298. https://doi.org/10.1145/1141885.1141892

[3]

Ernst Hairer e Gerhard Wanner. 1996. Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems (2ª ed.). Springer-Verlag, Berlin, Heidelberg. https://doi.org/10.1007/978-3-642-05221-7

---

1. Carl David Tolmé Runge (1856-1927) e Martin Wilhelm Kutta (1867-1944)

2. Carl Gottfried Neumann (1832-1925)

3. Karl Heun (1859-1929)

4. Anthony Ralston (1930-2019)

5. Brook Taylor (1685-1731)

6. Leonhard Euler (1707-1783)

7. Germund Dahlquist (1925-2005)

8. Carl Friedrich Gauss (1777-1855) e Adrien-Marie Legendre (1752-1833)

9. Jean Charles Rodolphe Radau (1835-1911)