Stabilità assoluta, problemi stiff e sistemi di equazioni differenziali

La nota precedente (Metodi di Eulero esplicito, implicito e Crank-NicolsonProblema di Cauchy $y'=f(t,y)$, $y(t_0)=y_0$, con $f$ lipschitziana in $y$ (costante $L$): ben posto. Nodi $t_n=t_0+nh$, $u_n\approx y(t_n)$. Eulero esplicito $u_{n+1}=u_n+h,f(t_n,u_n)$; Eulero implicito $u_{n+1}=u_n+h,f(t_{n+1},u_{n+1})$ (equazione da risolvere a ogni passo: punto fisso se $hL<1$, Newton, o formula chiusa se $f$ è lineare); Crank-Nicolson $u_{n+1}=u_n+\frac h2[f(t_n,u_n)+f(t_{n+1},u_{n+1})]$ (regola dei trapezi). Errore di troncamento locale $\tau_n=\frac{y(t_{n+1})-y(t_n)}h-\Phi$: $\frac h2y''$ (EE), $-\frac h2y''$ (EI), $-\frac{h^2}{12}y'''$ (CN). Consistenza + stabilità $\Rightarrow$ convergenza: ordine 1 per EE e EI, ordine 2 per CN, cioè errore globale $O(h)$ e $O(h^2)$; per EE $|e_n|\le\frac{Mh}{2L}\big(e^{L(t_n-t_0)}-1\big)$ con $M=\max|y''|$ (D). Dimezzando $h$ l'errore si divide per 2 (EE, EI) o per 4 (CN).Metodi di Eulero esplicito, implicito e Crank-Nicolson →) mostra che l'errore tende a zero per $h\to0$. Ma in pratica si usa un $h$ fisso, e il calcolatore ha un tempo finito. Un metodo di ordine alto con un passo troppo grande può dare risultati che esplodono in modo del tutto indipendente dall'accuratezza: è un problema di stabilità. Qui si studia quanto deve essere piccolo $h$ perché la soluzione numerica si comporti come quella vera (che decade, se il sistema fisico dissipa), cosa sono i problemi stiff, e come si trattano i sistemi di equazioni, che includono le equazioni di ordine superiore (l'oscillatore, il pendolo, i circuiti con più componenti).

Prerequisiti: Problema di Cauchy per sistemi di equazioni differenzialiUn sistema in forma normale è y' = f(t, y) con y a valori in Rⁿ; il problema di Cauchy aggiunge y(t₀) = y₀. Ogni equazione di ordine n si riduce a un sistema del primo ordine. Peano: f continua ⇒ esiste almeno una soluzione locale (non unica in generale: y' = y^(2/3)). Cauchy-Lipschitz: f e ∂f/∂y continue (o lipschitziana in y) ⇒ esistenza e unicità locali, e una soluzione massimale. Può esplodere in tempo finito (y' = y²); se |f(t,y)| ≤ H + K|y| su una striscia [a,b] × Rⁿ la soluzione è globale.Problema di Cauchy per sistemi di equazioni differenziali →, Sistemi lineari ed equazioni del secondo ordine a coefficienti costantiUn sistema lineare y' = A(t)y + g(t) con coefficienti continui ha, per ogni dato iniziale, una e una sola soluzione, definita su tutto l'intervallo (D). Le soluzioni del sistema omogeneo formano uno spazio vettoriale di dimensione n (D); quelle del sistema completo sono una soluzione particolare più l'integrale omogeneo (D). Per y'' + by' + cy = 0 si cercano e^(λt): tre casi (D). Un sistema 2×2 a coefficienti costanti si riduce a u'' − (tr A)u' + (det A)u = g, risolta con il metodo di somiglianza.Sistemi lineari ed equazioni del secondo ordine a coefficienti costanti →, Oscillatore massa-molla-smorzatore e transitoriMassa m, molla k e smorzatore viscoso c con forza esterna F(t): m x'' + c x' + k x = F(t), cioè x'' + (c/m) x' + (k/m) x = F/m. Grandezze: ω0 = √(k/m), smorzamento critico c_cr = 2√(km), ζ = c/c_cr. Libero: c > c_cr due esponenziali negativi (sovrasmorzato, nessuna oscillazione), c = c_cr (c1 + c2 t) e^(−ω0 t), c < c_cr e^(−ct/2m)(c1 cos ωd t + c2 sin ωd t) con ωd = ω0√(1−ζ²); con c > 0 ogni moto libero tende a 0 (transitorio). Forzato F0 cos Ωt con c > 0: x = transitorio + regime permanente X cos(Ωt − φ), X = F0/√((k−mΩ²)² + c²Ω²), tan φ = cΩ/(k−mΩ²); massimo di X per Ω² = k/m − c²/(2m²). Primo ordine m v' + c v = F: v = F/c + (v0 − F/c) e^(−ct/m), costante di tempo τ = m/c.Oscillatore massa-molla-smorzatore e transitori →, Autovalori e autovettoriUn autovettore è un vettore non nullo che una funzione lineare manda in un suo multiplo; si trovano gli autovalori come radici del polinomio caratteristico det(A − λI) e gli autovettori come nucleo di A − λI. Matrici simili hanno gli stessi autovalori.Autovalori e autovettori →, Fattorizzazione LU, pivoting e CholeskyRisolvere $A\mathbf x=\mathbf b$ con Gauss equivale a fattorizzare $A=LU$ ($L$ triangolare inferiore con $1$ sulla diagonale e i moltiplicatori $l_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$ sotto la diagonale, $U$ la matrice a scala finale) e poi risolvere $L\mathbf y=\mathbf b$ (in avanti) e $U\mathbf x=\mathbf y$ (all'indietro). Esiste (ed è unica) se i minori principali di testa $\det A_k$, $k<n$, sono $\ne0$. Costo $\frac23n^3$ flop per la fattorizzazione, $2n^2$ per ogni termine noto; $\det A=\prod u_{kk}$. Con un pivot piccolo l'algoritmo è instabile: il pivoting parziale (scambio di righe che porta in posizione di pivot l'elemento di modulo massimo della colonna) dà $PA=LU$ con $|l_{ik}|\le1$ ed esiste per ogni $A$ invertibile. Se $A$ è simmetrica definita positiva esiste $A=LL^T$ (Cholesky), con $l_{jj}=\sqrt{a_{jj}-\sum_{k<j}l_{jk}^2}$ e $l_{ij}=\frac{a_{ij}-\sum_{k<j}l_{ik}l_{jk}}{l_{jj}}$; costo $\frac13n^3$, senza pivoting.Fattorizzazione LU, pivoting e Cholesky → (i metodi impliciti per sistemi richiedono la risoluzione di un sistema lineare a ogni passo).

1. Il problema test e la stabilità assoluta

La soluzione di $y'=\lambda y$, $y(0)=1$, con $\lambda\in\mathbb C$ e $\mathrm{Re}\,\lambda<0$ è $e^{\lambda t}\to0$. Un metodo è assolutamente stabile con il passo $h$ se anche la soluzione numerica tende a $0$. Applicando i tre metodi (con $z=h\lambda$):

• EE: $u_{n+1}=u_n+h\lambda u_n=(1+z)\,u_n$;
• EI: $u_{n+1}=u_n+h\lambda u_{n+1}\Rightarrow u_{n+1}=\dfrac1{1-z}\,u_n$;
• CN: $u_{n+1}=u_n+\frac{h\lambda}2(u_n+u_{n+1})\Rightarrow u_{n+1}=\dfrac{1+z/2}{1-z/2}\,u_n$.

In tutti i casi $u_n=R(z)^n\,u_0$, quindi $u_n\to0$ se e solo se $|R(z)|<1$.

Definizione. La regione di stabilità assoluta è $\mathcal A=\{z\in\mathbb C:|R(z)|<1\}$. Il metodo è A-stabile se $\mathcal A$ contiene tutto il semipiano $\mathrm{Re}\,z<0$.

• EE. $|1+z|<1$: disco di centro $-1$ e raggio $1$. Per $\lambda<0$ reale: $|1-h|\lambda||<1\iff h<\dfrac2{|\lambda|}$. Condizionatamente stabile.
• EE con $\lambda=\mu+i\omega$ complesso ($\mu<0$): $|1+h\lambda|^2=(1+h\mu)^2+h^2\omega^2<1\iff 2h\mu+h^2(\mu^2+\omega^2)<0$, e dividendo per $h>0$ si ottiene $h<\dfrac{-2\mu}{\mu^2+\omega^2}=\dfrac{2\,|\mathrm{Re}\,\lambda|}{|\lambda|^2}$. Per $\omega=0$ si ritrova $\frac2{|\lambda|}$; con la parte immaginaria il limite è più severo, e per $\mu\to0$ (autovalori quasi immaginari, poco smorzamento) tende a zero. Esempio: $\lambda=-1\pm3i$ dà $h<\frac{2}{10}=0{,}2$, mentre $\frac2{|\lambda|}=0{,}63$ sarebbe sbagliato.
• EI. $|1-z|>1$: esterno del disco di centro $+1$ e raggio $1$; contiene tutto il semipiano sinistro (se $\mathrm{Re}\,z<0$ è $|1-z|>1$). A-stabile.
• CN. $|1+z/2|<|1-z/2|\iff\mathrm{Re}\,z<0$: la regione è esattamente il semipiano sinistro. A-stabile.

Perché conta. Se $\mathrm{Re}\,\lambda<0$ il sistema fisico è dissipativo e la soluzione decade. Con EE e $h$ troppo grande si ottiene una soluzione numerica che cresce in modo oscillante, fisicamente assurda. Con EI e CN non c'è nessuna restrizione di stabilità: $h$ si sceglie solo in base all'accuratezza.

Esempio ($\lambda=-50$, $y(0)=1$). Il fattore di EE è $1-50h$: $h=0{,}01$: $0{,}5$ (decade monotonamente); $h=0{,}02$: $0$ (azzera la soluzione al primo passo); $h=0{,}03$: $-0{,}5$ (oscilla smorzata); $h=0{,}04$: $-1$ (oscillazione permanente di ampiezza $1$); $h=0{,}05$: $-1{,}5$ (diverge). La soglia è $h=\frac{2}{50}=0{,}04$. Con EI il fattore è $\frac1{1+50h}$ ($h=0{,}05$: $0{,}286$) e con CN $\frac{1-25h}{1+25h}$ ($h=0{,}05$: $-0{,}111$): stabili.

Comportamento per $z\to-\infty$. $R_{EI}(z)=\frac1{1-z}\to0$: l'EI smorza immediatamente le componenti molto veloci (come fa la soluzione vera). $R_{CN}(z)\to-1$: la CN le fa decadere molto lentamente e con segno alternato: per $\lambda=-100$, $h=0{,}1$ ($z=-10$): $R_{CN}=\frac{1-5}{1+5}=-0{,}667$, $R_{EI}=\frac1{11}=0{,}091$, $R_{EE}=-9$. La CN è accurata ma può produrre oscillazioni residue nei problemi stiff con passi grandi.

2. Problemi stiff

Un problema è stiff (rigido) quando contiene scale temporali molto diverse: una parte della soluzione che varia rapidamente (un transitorio che muore in fretta) e una che varia lentamente.

Esempio. $y'=-1000\,(y-\cos t)-\sin t$, $y(0)=0$. La soluzione è $y(t)=\cos t-e^{-1000t}$: dopo $\approx0{,}005$ s il transitorio $e^{-1000t}$ è sparito e $y\approx\cos t$ è liscia, adatta a passi di $0{,}1$ s. Ma la stabilità di Eulero esplicito richiede $h<\frac2{1000}=0{,}002$ per tutta la durata, anche quando il transitorio non c'è più: per 10 s sono $5000$ passi invece dei $100$ che l'accuratezza richiederebbe. Con EI o CN si usa $h=0{,}1$.

Per un sistema $\mathbf y'=A\mathbf y$ lo stiffness ratio è $\dfrac{\max|\mathrm{Re}\,\lambda_i|}{\min|\mathrm{Re}\,\lambda_i|}$. Esempi in ingegneria: reti di reazioni chimiche con costanti molto diverse, circuiti con $RC$ molto diversi, un sistema meccanico con molla molto rigida e una massa pesante, il raffreddamento di un pezzo che ha una superficie sottile e un nucleo massiccio.

Rimedi. I metodi impliciti (EI, CN e quelli di ordine più alto) non hanno limiti di stabilità. Il costo è il sistema lineare (o non lineare, con Newton) da risolvere a ogni passo. In MATLAB, ode45 (esplicito) con un problema stiff usa passi minuscoli e va lentissimo; si cambia in ode15s o ode23s. Un segnale pratico di stiffness è proprio questo: il solutore esplicito fa passi molto più piccoli di quanto la soluzione richieda.

3. Sistemi di equazioni e equazioni di ordine superiore

Un sistema $\mathbf y'=\mathbf f(t,\mathbf y)$, $\mathbf y(t_0)=\mathbf y_0$ si risolve componente per componente con le stesse formule, trattando $\mathbf y$ come vettore: $$\text{EE: }\mathbf u_{n+1}=\mathbf u_n+h\,\mathbf f(t_n,\mathbf u_n),\qquad\text{EI: }\mathbf u_{n+1}=\mathbf u_n+h\,\mathbf f(t_{n+1},\mathbf u_{n+1}),\qquad\text{CN: }\mathbf u_{n+1}=\mathbf u_n+\tfrac h2\big[\mathbf f_n+\mathbf f_{n+1}\big].$$ Consistenza, ordine e convergenza valgono con le norme di vettore (la costante di Lipschitz è calcolata con una norma).

Riduzione all'ordine uno. Un'equazione $y^{(m)}=F(t,y,y',\dots,y^{(m-1)})$ equivale al sistema per $\mathbf y=(y,y',\dots,y^{(m-1)})$. L'oscillatore smorzato $m\,x''+c\,x'+k\,x=F(t)$ con $v=x'$: $$\begin{pmatrix}x\\v\end{pmatrix}'=\begin{pmatrix}0&1\\-\frac km&-\frac cm\end{pmatrix}\begin{pmatrix}x\\v\end{pmatrix}+\begin{pmatrix}0\\\frac{F(t)}m\end{pmatrix}.$$ Gli autovalori sono $\lambda=\dfrac{-c\pm\sqrt{c^2-4mk}}{2m}$.

Sistema lineare $\mathbf y'=A\mathbf y$.

• EE: $\mathbf u_{n+1}=(I+hA)\,\mathbf u_n$.
• EI: $(I-hA)\,\mathbf u_{n+1}=\mathbf u_n$: si risolve un sistema lineare; la matrice $I-hA$ è sempre la stessa se $h$ è costante, quindi si fa una fattorizzazione LU all'inizio e a ogni passo solo due sostituzioni ($2m^2$ operazioni).
• CN: $\big(I-\tfrac h2A\big)\,\mathbf u_{n+1}=\big(I+\tfrac h2A\big)\,\mathbf u_n$.

Stabilità per sistemi. Se $A$ è diagonalizzabile, $A=V\Lambda V^{-1}$, il sistema si disaccoppia in $m$ equazioni test $y_i'=\lambda_iy_i$ nelle coordinate $V^{-1}\mathbf y$. Quindi:

Stabilità per sistemi. Il metodo è assolutamente stabile se $h\lambda_i\in\mathcal A$ per ogni autovalore $\lambda_i$ di $A$.

Per EE con autovalori reali negativi serve $h<\dfrac2{\max_i|\lambda_i|}$: l'autovalore più grande in modulo (cioè la componente più veloce) fissa il passo.

Esempio (sistema stiff). $A=\begin{pmatrix}-50{,}5&49{,}5\\49{,}5&-50{,}5\end{pmatrix}$ ha autovalori $-1$ (autovettore $(1,1)$) e $-100$ (autovettore $(1,-1)$). Con $\mathbf y(0)=(1,0)$ la soluzione è $\mathbf y=\frac12e^{-t}(1,1)+\frac12e^{-100t}(1,-1)$: già per $t=0{,}05$ resta solo $\frac12e^{-t}(1,1)$. Stabilità di EE: $h<\frac2{100}=0{,}02$. Valore di $y_1(1)$ (esatto $0{,}18394$):

Per $h=0{,}02$ EE è al limite ($|1-100h|=1$): la componente veloce non decade e il risultato è completamente sbagliato. Per $h\ge0{,}025$ esplode. EI è stabile e perde precisione solo gradualmente (ordine 1); CN è molto accurata finché $|h\lambda|$ non è troppo grande: a $h=0{,}1$ ($z=-10$) il fattore $-0{,}667$ della componente veloce non è ancora morto e dà un errore $\approx0{,}009$.

4. Oscillatore: ampiezza e energia

Per un moto senza smorzamento, $x''+\omega^2x=0$, gli autovalori sono immaginari puri $\lambda=\pm i\omega$, sull'asse immaginario (bordo del semipiano). La soluzione esatta conserva l'energia $E=\frac12(v^2+\omega^2x^2)$. Con $z=\pm i\omega h$:

• EE: $|R|=|1+i\omega h|=\sqrt{1+\omega^2h^2}>1$ per qualunque $h$. L'energia (proporzionale al quadrato dell'ampiezza) cresce di un fattore $1+\omega^2h^2$ a ogni passo: instabile sempre. Nessun $h$ lo rende accettabile su tempi lunghi (si può solo ridurre il tasso di crescita).
• EI: $|R|=\frac1{|1-i\omega h|}=\frac1{\sqrt{1+\omega^2h^2}}<1$: l'energia decade di un fattore $\frac1{1+\omega^2h^2}$ per passo (smorzamento numerico: l'oscillazione sembra smorzata anche se non lo è).
• CN: $|R|=\frac{|1+i\omega h/2|}{|1-i\omega h/2|}=1$ esattamente: l'energia è conservata, c'è solo un piccolo errore di fase.

Esempio ($\omega=1$, $x(0)=1$, $v(0)=0$, $h=0{,}1$). Rapporto $E_n/E_0$:

(EE: $(1+h^2)^n=1{,}01^n$; EI: $1{,}01^{-n}$.) Dopo un periodo ($t\approx6{,}3$) l'energia con EE è quasi raddoppiata e con EI è dimezzata. In ingegneria questo conta: simulando un sistema meccanico oscillante per molti periodi, Eulero esplicito aggiunge energia (e un sistema può «esplodere»), Eulero implicito ne toglie (e sembra frenato). Per la CN l'errore è solo nella frequenza: l'argomento di $R$ è $2\arctan\frac{\omega h}2\approx\omega h-\frac{(\omega h)^3}{12}$ contro $\omega h$ del moto esatto, quindi la fase si sposta con errore $O(h^2)$.

5. Implementazione

import numpy as np

def integra(A, y0, h, n, metodo):
    """y' = A y con Eulero esplicito (ee), implicito (ei) o Crank-Nicolson (cn). Restituisce y_n."""
    I = np.eye(len(y0))
    if metodo == "ee":
        P = I + h * A                                      # y_{n+1} = (I + hA) y_n
    elif metodo == "ei":
        P = np.linalg.inv(I - h * A)                       # (I - hA) y_{n+1} = y_n
    else:
        P = np.linalg.solve(I - h / 2 * A, I + h / 2 * A)  # (I - hA/2) y_{n+1} = (I + hA/2) y_n
    return np.linalg.matrix_power(P, n) @ y0

# 1) oscillatore x'' + x = 0 -> y = (x, v), y' = A y, energia E = (x^2 + v^2)/2
A = np.array([[0.0, 1.0], [-1.0, 0.0]])
y0 = np.array([1.0, 0.0])
print("oscillatore, h = 0.1, energia E/E0 dopo n passi:")
for n in (10, 63, 200):
    E = [np.sum(integra(A, y0, 0.1, n, m) ** 2) / np.sum(y0 ** 2) for m in ("ee", "ei", "cn")]
    print(f"  n = {n:3d} (t = {0.1 * n:5.1f}):  EE {E[0]:.4f}   EI {E[1]:.4f}   CN {E[2]:.4f}")

# 2) sistema stiff con autovalori -1 e -100
A = np.array([[-50.5, 49.5], [49.5, -50.5]])
y0 = np.array([1.0, 0.0])
esatta = 0.5 * np.exp(-1.0) * np.array([1, 1]) + 0.5 * np.exp(-100.0) * np.array([1, -1])
print("\nsistema stiff, y(1) esatta =", esatta.round(6))
for h in (0.01, 0.02, 0.025, 0.05, 0.1):
    n = round(1 / h)
    print(f"  h = {h:5.3f}:", *[f"{m.upper()} {integra(A, y0, h, n, m)[0]: .4e}" for m in ("ee", "ei", "cn")])

Output:

oscillatore, h = 0.1, energia E/E0 dopo n passi:
  n =  10 (t =   1.0):  EE 1.1046   EI 0.9053   CN 1.0000
  n =  63 (t =   6.3):  EE 1.8717   EI 0.5343   CN 1.0000
  n = 200 (t =  20.0):  EE 7.3160   EI 0.1367   CN 1.0000

sistema stiff, y(1) esatta = [0.18394 0.18394]
  h = 0.010: EE  1.8302e-01 EI  1.8486e-01 CN  1.8394e-01
  h = 0.020: EE  6.8208e-01 EI  1.8576e-01 CN  1.8393e-01
  h = 0.025: EE  5.5287e+06 EI  1.8622e-01 CN  1.8393e-01
  h = 0.050: EE  5.4976e+11 EI  1.8844e-01 CN  1.8390e-01
  h = 0.100: EE  1.7434e+09 EI  1.9277e-01 CN  1.9246e-01

(La potenza matrix_power equivale a ripetere il passo $n$ volte; in un programma vero si fa un ciclo e si memorizza ogni passo.) Per sistemi non lineari si usa Newton a ogni passo per gli schemi impliciti. Equivalenti MATLAB: ode45 (non stiff, esplicito Runge-Kutta), ode15s, ode23s (stiff), [t,y] = ode15s(@(t,y) A*y, [0 1], [1;0]); i metodi di questa nota si programmano con un ciclo e (eye(2)-h*A)\y.

Errori da evitare

• Scambiare accuratezza e stabilità: un ordine alto non dà stabilità; Eulero esplicito è consistente di ordine 1 ma instabile per $h>\frac2{|\lambda|}$.
• Dimenticare di controllare TUTTI gli autovalori: basta uno fuori dalla regione perché la soluzione numerica esploda, anche se la sua componente nella soluzione vera è trascurabile.
• Usare $h<\frac2{|\lambda|}$ per autovalori complessi: il limite $\frac2{|\lambda|}$ vale solo per $\lambda$ reale negativo; per $\lambda=\pm i\omega$ EE è instabile per ogni $h$, e per $\lambda=\mu\pm i\omega$ la condizione è $|1+h\lambda|<1$ con modulo complesso, cioè $h<\frac{2|\mu|}{\mu^2+\omega^2}$.
• Credere che CN risolva tutto nei problemi stiff: per $z\to-\infty$ è $R\to-1$, quindi con $h$ molto grande restano oscillazioni residue; per problemi molto stiff si preferiscono metodi impliciti con $R(\infty)=0$ (come EI e gli schemi BDF).
• Calcolare l'inversa $(I-hA)^{-1}$ a ogni passo: si fattorizza una volta e si risolve.
• Pensare che lo smorzamento numerico di EI sia innocuo: sopprime oscillazioni reali (energia dissipata che non c'è).

Esercizi svolti: Esercizio - passo massimo di Eulero esplicito per un'equazione test e confronto con i metodi impliciti, Esercizio - sistema 2x2 stiff con autovalori -1 e -100 e limite di stabilità, Esercizio - oscillatore armonico con Eulero esplicito, implicito e Crank-Nicolson, energia, Esercizio - massa-molla-smorzatore ridotto al primo ordine e passo di stabilità.

Domande d'esame

1. Definire la regione di stabilità assoluta e determinarla per Eulero esplicito, implicito e Crank-Nicolson. Traccia. Problema test $y'=\lambda y$, $u_{n+1}=R(z)u_n$; $|1+z|<1$; $|1-z|>1$; $\mathrm{Re}\,z<0$.
2. Che cosa significa A-stabile? Quali dei tre metodi lo sono? Perché EE non può esserlo? Traccia. $\mathcal A\supset\{\mathrm{Re}\,z<0\}$; EI e CN; per EE $\mathcal A$ è un disco limitato.
3. Che cos'è un problema stiff? Perché Eulero esplicito è inadatto? Traccia. Scale di tempo molto diverse; $h<\frac2{|\lambda|_{\max}}$ imposto dalla componente veloce anche quando è già estinta.
4. Per l'oscillatore $x''+\omega^2x=0$ studiare come variano ampiezza ed energia per i tre metodi. Traccia. $|R_{EE}|=\sqrt{1+\omega^2h^2}>1$, $|R_{EI}|<1$, $|R_{CN}|=1$.