Condizionamento di un problema e stabilità di un algoritmo

Prerequisiti: Aritmetica di macchina - floating point, precisione di macchina e cancellazioneUn calcolatore rappresenta solo un insieme finito di numeri F(β,t,L,U); ogni numero reale viene arrotondato con errore relativo al più u = ½β^(1−t) (u ≈ 1,1·10⁻¹⁶ in doppia precisione) e ogni operazione è esatta più un errore (1+δ), |δ| ≤ u. Conseguenze: overflow e underflow, somma non associativa, assorbimento, e soprattutto la cancellazione numerica nella differenza di numeri quasi uguali, che si evita riformulando l'espressione.Aritmetica di macchina - floating point, precisione di macchina e cancellazione → (errore relativo, unità di arrotondamento $u$), Derivata - definizione e significatoLa derivata è il limite del rapporto incrementale; geometricamente è la pendenza della retta tangente. f è derivabile in x0 se e solo se f(x) = f(x0) + f'(x0)(x − x0) + o(x − x0); derivabile implica continua, non viceversa. Derivata destra e sinistra, punti angolosi, flessi a tangente verticale, cuspidi.Derivata - definizione e significato → e Teoremi di Rolle e di LagrangeRolle: f continua su [a,b], derivabile in ]a,b[, f(a) = f(b) ⇒ esiste c in ]a,b[ con f'(c) = 0. Lagrange (valor medio): stesse ipotesi senza f(a) = f(b) ⇒ esiste c con f'(c) = (f(b) − f(a))/(b − a): la tangente in c è parallela alla corda. Conseguenza: su un intervallo f è costante se e solo se f' = 0; su un insieme che non è un intervallo è solo costante a tratti (arctan x + arctan(1/x)).Teoremi di Rolle e di Lagrange →.

Il calcolo numerico si occupa di due domande distinte che è facile confondere:

• quanto è sensibile il problema ai piccoli errori sui dati? (condizionamento, proprietà del problema);
• quanto è sensibile il procedimento con cui lo risolviamo agli errori di arrotondamento? (stabilità, proprietà dell'algoritmo).

Le sorgenti dell'errore

L'errore con cui si ottiene un risultato ha tre origini:

1. Errore sui dati (o inerente): i dati di partenza sono misure (un calibro dà $12{,}347\,\text{mm}$ con incertezza $\pm 5\cdot 10^{-4}\,\text{mm}$) oppure numeri reali arrotondati per entrare nel calcolatore.
2. Errore di troncamento (o analitico): il metodo sostituisce un oggetto infinito o continuo con uno finito. Esempio: si calcola una derivata con un rapporto incrementale, si arresta una serie, si approssima un integrale con una somma.
3. Errore di arrotondamento: ogni operazione di macchina commette un errore relativo $\le u$.

Problemi ben posti

Definizione (problema ben posto, Hadamard). Il problema "dato $d$, trova $x$ tale che $F(x, d) = 0$" è ben posto se: (1) la soluzione $x$ esiste; (2) è unica; (3) dipende con continuità dai dati $d$. Altrimenti è mal posto.

La terza condizione è la più importante in pratica: se una piccola variazione di $d$ può produrre una grande variazione di $x$, i dati inevitabilmente sporchi (misure, arrotondamenti) rendono il risultato inaffidabile. La continuità però non basta: serve sapere quanto il problema amplifica.

Numero di condizionamento

Si perturba il dato $x$ in $\tilde x = x(1 + \varepsilon)$, con $|\varepsilon|$ piccolo, e si calcola $y = f(x)$. Per la formula di Lagrange (o di Taylor al primo ordine, Formula di Taylor con resto di LagrangeSe f ha n+1 derivate continue su [a,b], l'errore tra f e il suo polinomio di Taylor si scrive esattamente come f^(n+1)(c)(x-x0)^(n+1)/(n+1)! per un c opportuno: serve a stimare l'errore su tutto un intervallo e a scegliere il grado del polinomio.Formula di Taylor con resto di Lagrange →) $$f(\tilde x) - f(x) = f'(\xi)(\tilde x - x) \approx f'(x)\,x\,\varepsilon,$$ e dividendo per $f(x)$ si ottiene l'errore relativo sul risultato: $$\frac{f(\tilde x) - f(x)}{f(x)} \approx \underbrace{\frac{x f'(x)}{f(x)}}_{K(x)}\;\varepsilon.$$

Definizione (numero di condizionamento). Per $y = f(x)$ derivabile con $x f(x) \ne 0$ $$K(x) = \left|\frac{x\,f'(x)}{f(x)}\right|$$ è il fattore con cui un errore relativo sul dato $x$ viene moltiplicato nel risultato. Il problema è ben condizionato se $K$ è piccolo (qualche unità), mal condizionato se $K$ è grande.

Regola pratica. Se i dati hanno $k$ cifre significative corrette, il risultato ne ha circa $k - \log_{10} K$. Un errore di solo arrotondamento ($u \approx 10^{-16}$) con $K = 10^{10}$ lascia solo $6$ cifre.

Altri casi con più variabili: per $y = f(x_1, \dots, x_n)$ ogni dato ha il suo $K_i = \left|\frac{x_i\,\partial_i f}{f}\right|$ e l'errore massimo è circa $\sum_i K_i |\varepsilon_i|$.

Esempio 1: il coseno vicino a $\frac\pi2$. Per $x = 1{,}57$: $K = 1{,}57 \cdot \tan 1{,}57 = 1{,}57 \cdot 1255{,}77 \approx 1971{,}6$. Si perturba $x$ con errore relativo $10^{-4}$, cioè $\tilde x = 1{,}57 \cdot 1{,}0001 = 1{,}570157$. Previsione: errore relativo su $\cos x$ circa $1971{,}6 \cdot 10^{-4} \approx 0{,}197$ ($19{,}7\%$). Calcolo diretto: $\cos 1{,}57 = 7{,}9633\cdot 10^{-4}$ e $\cos 1{,}570157 = 6{,}3933\cdot 10^{-4}$, variazione relativa $\frac{1{,}570\cdot 10^{-4}}{7{,}9633\cdot 10^{-4}} = 0{,}1972$ ✓. Una perturbazione dello $0{,}01\%$ sull'angolo cambia il risultato del $20\%$: vicino a $\frac\pi2$ la pendenza del coseno vale circa $-1$, quindi la variazione assoluta del risultato è pari a quella del dato ($1{,}57\cdot 10^{-4}$); ma il valore del coseno è prossimo a zero, e diviso per un numero così piccolo l'errore relativo esplode.

Esempio 2: radice doppia. Le radici di $(x-1)^2 = \varepsilon$ sono $x = 1 \pm \sqrt\varepsilon$. Se il termine noto del polinomio $x^2 - 2x + 1$ viene perturbato di $\varepsilon = 10^{-10}$ (errore minuscolo, come un arrotondamento) le radici si spostano di $\sqrt\varepsilon = 10^{-5}$: amplificazione $10^{5}$. Per $\varepsilon = -10^{-10}$ le radici diventano complesse $1 \pm 10^{-5}i$ (Python: np.roots([1, -2, 1 + 1e-10]) dà 1 ± 1e-5j). Una radice multipla è sempre mal condizionata: il grafico del polinomio è tangente all'asse, una piccola traslazione verticale lo allontana molto in orizzontale.

Esempio 3: sistema lineare. Il sistema $x + y = 2$, $x + 1{,}0001\,y = 2{,}0001$ ha soluzione $(1, 1)$. Se il secondo termine noto passa da $2{,}0001$ a $2{,}0002$ (variazione relativa $5\cdot 10^{-5}$) la soluzione diventa $(0, 2)$: il risultato cambia del $100\%$. Le due rette sono quasi parallele; amplificazione dell'ordine di $2\cdot 10^4$. Il numero di condizionamento di un sistema è trattato in Norme di matrici e condizionamento dei sistemi lineariNorme di vettore $|x|1=\sum|x_i|$, $|x|2=\sqrt{\sum x_i^2}$, $|x|\infty=\max|x_i|$. Norma di matrice indotta $|A|=\max{x\ne0}\frac{|Ax|}{|x|}$: $|A|1$ = massima somma per colonne, $|A|\infty$ = massima somma per righe, $|A|2=\sqrt{\lambda{\max}(A^TA)}$ (= $\rho(A)$ se $A$ simmetrica); valgono $|Ax|\le|A||x|$, $|AB|\le|A||B|$, $\rho(A)\le|A|$. Numero di condizionamento $\kappa(A)=|A||A^{-1}|\ge1$ ($\kappa_2=\frac{\sigma_{\max}}{\sigma_{\min}}$, per SPD $\frac{\lambda_{\max}}{\lambda_{\min}}$, per ortogonali 1). Se $A\mathbf x=\mathbf b$ e $A(\mathbf x+\delta\mathbf x)=\mathbf b+\delta\mathbf b$: $\frac{|\delta\mathbf x|}{|\mathbf x|}\le\kappa(A)\frac{|\delta\mathbf b|}{|\mathbf b|}$ (D); con residuo $\mathbf r=\mathbf b-A\hat{\mathbf x}$: $\frac1\kappa\frac{|\mathbf r|}{|\mathbf b|}\le\frac{|\mathbf x-\hat{\mathbf x}|}{|\mathbf x|}\le\kappa\frac{|\mathbf r|}{|\mathbf b|}$. Si perdono circa $\log_{10}\kappa$ cifre; un residuo piccolo non garantisce un errore piccolo se $\kappa$ è grande.Norme di matrici e condizionamento dei sistemi lineari →.

Stabilità di un algoritmo

Un algoritmo calcola $y = f(x)$ con una successione di operazioni elementari; ciascuna introduce un errore $\le u$, e questi errori si propagano e si possono amplificare. Sia $\hat y$ il risultato calcolato.

Definizione (errore in avanti e all'indietro).

• Errore in avanti (forward): $\dfrac{|\hat y - f(x)|}{|f(x)|}$.
• Errore all'indietro (backward): la perturbazione $\Delta x$ dei dati per cui $\hat y$ è il risultato esatto, $\hat y = f(x + \Delta x)$, cioè $\dfrac{|\Delta x|}{|x|}$. Un algoritmo è stabile all'indietro se l'errore all'indietro è dell'ordine di $u$ (di poco più di un arrotondamento sui dati).

Collegamento fondamentale, dalla definizione di $K$: $$\text{errore in avanti} \;\lesssim\; K(x) \cdot \text{errore all'indietro} \;\approx\; K(x)\cdot u .$$ Quindi un algoritmo stabile applicato a un problema ben condizionato dà un risultato accurato. Se invece il problema è mal condizionato ($K$ grande), nessun algoritmo può dare di più di circa $\log_{10}(1/(K u))$ cifre corrette: l'errore è nel problema, non nel metodo. Se il problema è ben condizionato e il risultato è sbagliato, l'algoritmo è instabile.

Esempio svolto: ricorrenza instabile. Sia $I_n = \displaystyle\int_0^1 x^n e^{x-1}\,dx$ per $n \ge 0$. Per parti (Integrazione per parti∫ f·g' dx = f·g − ∫ f'·g dx (regola del prodotto letta al contrario). Si usa per log x, arcsin x, arctan x (scritti come "funzione per 1"), per sin²x, cos²x, sinh²x, cosh²x (integrale circolare: l'integrale di partenza ricompare e si porta a sinistra) e per prodotti come x^n·e^(αx), x^n·sin(βx), e^(αx)·sin(βx).Integrazione per parti →): $$I_n = \bigl[x^n e^{x-1}\bigr]_0^1 - n\int_0^1 x^{n-1}e^{x-1}dx = 1 - n\,I_{n-1}, \qquad I_0 = 1 - e^{-1} = 0{,}63212\ldots$$ Il problema è ben condizionato (un integrale di una funzione positiva e limitata: $\frac{1}{e(n+1)} \le I_n \le \frac1{n+1}$, cioè $I_n$ è piccolo e decrescente). Si prova la ricorrenza in avanti con $4$ cifre significative, partendo da $I_0 \to 0{,}6321$ (errore $e_0 = 2{,}06\cdot 10^{-5}$):

Il motivo: se $\tilde I_{n-1} = I_{n-1} + e_{n-1}$, allora $\tilde I_n = 1 - n(I_{n-1} + e_{n-1})$ e quindi $e_n = -n\,e_{n-1}$, cioè $$|e_n| = n!\,|e_0|.$$ Ad ogni passo l'errore viene moltiplicato per $n$: già a $n = 7$ il risultato è sbagliato quasi del $100\%$. In doppia precisione $e_0 \approx 10^{-16}$, e $18! \cdot 10^{-16} \approx 0{,}6$: la ricorrenza in avanti dà $I_{20} = -30{,}19$, un valore assurdo (l'integrale di una funzione positiva è positivo).

Rimedio. Si inverte la ricorrenza: $I_{n-1} = \dfrac{1 - I_n}{n}$. Ora l'errore viene diviso per $n$ a ogni passo: $e_{n-1} = -e_n/n$. Si può partire da un valore anche grossolano a un indice alto $N$, ad esempio $I_N \approx 0$ (l'errore è $\le \frac{1}{N+1}$), e a $n$ basso l'errore residuo è $\frac{n!}{N!}\cdot\frac{1}{N+1}$. Con $4$ cifre e partendo da $I_{10} \approx 0$ si ottiene all'indietro $I_9 = 0{,}1000$ (esatto $0{,}0916$), $I_8 = 0{,}1000$ ($0{,}1009$), $I_7 = 0{,}1125$ ($0{,}1124$), $I_6 = 0{,}1268$ ($0{,}1268$), $\dots$, $I_0 = 0{,}6321$ ✓: le cifre corrette sono quelle del dato, nessuna amplificazione. Lo stesso problema, due algoritmi: uno instabile, uno stabile.

Consistenza, stabilità e convergenza di uno schema numerico

Un metodo numerico dipende da un parametro che ne fissa l'accuratezza (il passo $h$, il numero $n$ di iterazioni o nodi): $F_n(x_n, d) = 0$ approssima $F(x, d) = 0$.

Definizione.

• Consistenza: se si inserisce nello schema la soluzione esatta $x$, il residuo tende a zero per $h \to 0$ ($n \to \infty$): lo schema "assomiglia" al problema.
• Stabilità: la soluzione numerica dipende con continuità, in modo uniforme rispetto a $h$, dai dati: piccole perturbazioni dei dati non vengono amplificate al diminuire di $h$.
• Convergenza: $x_n \to x$ per $h \to 0$.

Teorema (di equivalenza di Lax-Richtmyer). Per un problema ben posto e lineare, uno schema consistente è convergente se e solo se è stabile. Per i problemi non lineari vale in generale almeno la direzione consistenza + stabilità $\Rightarrow$ convergenza.

Perché funziona. La soluzione esatta $x$, inserita nello schema, non lo soddisfa esattamente: lascia un residuo $\tau_h$ (errore di troncamento), cioè $x$ è la soluzione dello schema con dati perturbati di $\tau_h$. La soluzione numerica $x_h$ è invece la soluzione dello schema con i dati veri. La stabilità dice che due soluzioni dello schema con dati che differiscono di $\tau_h$ distano al più $C\,|\tau_h|$, con $C$ indipendente da $h$. Quindi $|x - x_h| \le C\,|\tau_h| \to 0$ per la consistenza. Se $C$ dipendesse da $h$ e crescesse (per esempio come $1/h$), il prodotto $C\,|\tau_h|$ potrebbe non tendere a zero: per questo la stabilità deve essere uniforme. Lo schema si applica alle equazioni differenziali in 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 →, dove $\tau_h$ è l'errore di troncamento locale.

Esempio svolto: derivata numerica. Per approssimare $f'(x)$ si usa il rapporto incrementale $D_h = \dfrac{f(x + h) - f(x)}{h}$.

• Errore di troncamento: dallo sviluppo di Taylor, $D_h - f'(x) = \frac h2 f''(\xi)$. Tende a zero con $h$ (schema consistente, di ordine $1$).
• Errore di arrotondamento: i due valori $f(x+h)$ e $f(x)$ hanno ciascuno errore relativo $u$; la loro differenza (numeri quasi uguali!) ha errore assoluto $\approx 2u|f|$, e dividendo per $h$ l'errore diventa $\approx \frac{2u|f|}{h}$: cresce al diminuire di $h$.
• Errore totale $E(h) \approx \dfrac h2|f''| + \dfrac{2u|f|}{h}$. Il minimo è per $h_{\text{ott}} = 2\sqrt{\dfrac{u|f|}{|f''|}}$, circa $\sqrt u \approx 10^{-8}$ in doppia precisione.

Per $f = e^x$ in $x = 1$ ($f' = e = 2{,}71828\ldots$), errore relativo calcolato:

L'errore scende come $h$ (troncamento) fino a $h \approx 10^{-8}$, poi risale come $1/h$ (arrotondamento): rimpicciolire $h$ oltre il necessario peggiora il risultato. Nel linguaggio delle definizioni: lo schema è consistente (il troncamento tende a zero), ma una perturbazione $\delta$ dei dati $f(x)$, $f(x+h)$ produce sul risultato un errore fino a $2\delta/h$, con un fattore di amplificazione $2/h$ che non è limitato uniformemente in $h$. Lo schema non è stabile rispetto agli errori di arrotondamento, e infatti in aritmetica di macchina non converge per $h \to 0$.

Codice

import math

I = 1 - math.exp(-1)                     # I_n = integrale di x^n e^(x-1) su [0,1]
for n in range(1, 21):
    I = 1 - n * I                        # in avanti: instabile
print("avanti   I_20 =", I)
J = 0.0
for n in range(40, 20, -1):
    J = (1 - J) / n                      # all'indietro: stabile, parte da I_40 ~ 0
print("indietro I_20 =", J)

for k in (2, 4, 6, 8, 10, 12):           # derivata numerica di exp in x = 1
    h = 10.0**-k
    print(k, abs((math.exp(1 + h) - math.exp(1)) / h - math.e) / math.e)

Output:

avanti   I_20 = -30.19239488558378
indietro I_20 = 0.045544884075818054
2 0.005016708416805288
4 5.000166739738369e-05
6 4.999377015549417e-07
8 2.429016271026634e-09
10 5.693704999536528e-07
12 0.00015903952213936482

In MATLAB il numero di condizionamento di una matrice è cond(A); per le derivate diff e per le radici di un polinomio roots.

Errori tipici

• Confondere condizionamento e stabilità: il primo riguarda il problema (indipendente dall'algoritmo), la seconda l'algoritmo. Un metodo instabile si cambia; un problema mal condizionato si può solo riformulare o accettare.
• Credere che ridurre $h$ migliori sempre un'approssimazione: oltre $h_{\text{ott}}$ domina l'arrotondamento.
• Calcolare $K$ come $|f'(x)|$ (condizionamento assoluto): il fattore relativo è $\left|\frac{xf'(x)}{f(x)}\right|$.
• Dedurre dal fatto che l'algoritmo è "corretto" (formula giusta) che sia anche stabile: la ricorrenza $I_n = 1 - nI_{n-1}$ è esatta ma instabile.

Domande d'esame

1. Quando un problema è ben posto? Risposta: se la soluzione esiste, è unica e dipende con continuità dai dati.
2. Definire il numero di condizionamento di $y = f(x)$ e calcolarlo per $\sqrt x$. Risposta: $K = |xf'(x)/f(x)|$; per $\sqrt x$ vale $x\cdot\frac{1}{2\sqrt x}/\sqrt x = \frac12$.
3. Differenza tra problema mal condizionato e algoritmo instabile. Risposta: il primo amplifica gli errori sui dati qualunque sia il metodo; il secondo amplifica gli errori di arrotondamento anche su un problema ben condizionato.
4. Perché nella derivata numerica esiste un $h$ ottimale? Risposta: il troncamento $\propto h$ cala con $h$, l'arrotondamento $\propto u/h$ cresce: la somma ha un minimo a $h \approx \sqrt u$.