Sistemi lineari al computer: LU, pivoting, Cholesky e metodi iterativi
Fattorizzazione LU e suo costo, pivoting parziale con esempio, Cholesky, matrici sparse, metodi di Jacobi e Gauss-Seidel, raggio spettrale e residuo.
In Geometria e Algebra i sistemi lineari si risolvevano con il metodo di Gauss su matrici 3×3. Nei problemi reali le matrici hanno migliaia o milioni di righe: la simulazione di un circuito integrato, l’analisi strutturale a elementi finiti, la ricostruzione di un’immagine in una TAC. Qui contano il costo del calcolo, la stabilità rispetto agli arrotondamenti e la struttura della matrice. Vediamo i metodi diretti (fattorizzazioni LU e di Cholesky) e quelli iterativi (Jacobi e Gauss-Seidel).
La fattorizzazione LU
L’eliminazione di Gauss, registrando i moltiplicatori usati, scompone la matrice nel prodotto
![\[ A = L\,U \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-e5ac2901a64b8e8d48bd11ca15341b44_l3.png "Rendered by QuickLaTeX.com")
con 
triangolare inferiore (con 1 sulla diagonale e i moltiplicatori sotto) e 
triangolare superiore (la forma a scala). Il sistema 
diventa due sistemi triangolari, che si risolvono per sostituzione:
![\[ L\mathbf y = \mathbf b \;\;(\text{in avanti}),\qquad U\mathbf x = \mathbf y \;\;(\text{all'indietro}) \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-e93d399d2e365b117570617fedad1f78_l3.png "Rendered by QuickLaTeX.com")
Il vantaggio. La fattorizzazione costa circa 
operazioni, ma le due sostituzioni solo 
. Se si devono risolvere molti sistemi con la stessa matrice e termini noti diversi (come in un simulatore che a ogni passo cambia solo le sorgenti), si fattorizza una volta sola. Con 
: 
operazioni per la fattorizzazione, 
per ogni nuovo termine noto.
Esempio.
![\[ A = \begin{pmatrix} 2 & 1 & 1\\ 4 & 3 & 3\\ 8 & 7 & 9\end{pmatrix} = \begin{pmatrix} 1 & 0 & 0\\ 2 & 1 & 0\\ 4 & 3 & 1\end{pmatrix}\begin{pmatrix} 2 & 1 & 1\\ 0 & 1 & 1\\ 0 & 0 & 2\end{pmatrix} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-84742d209aa5acdd0a8e0d47ed364837_l3.png "Rendered by QuickLaTeX.com")
I moltiplicatori 2 e 4 eliminano la prima colonna, il 3 la seconda.
Il pivoting
Se durante l’eliminazione compare uno zero sulla diagonale il metodo si blocca; se compare un numero molto piccolo, si divide per esso e gli errori di arrotondamento esplodono.
Esempio. Con 
e 
, la soluzione vera è circa 
. Senza scambi, il moltiplicatore è 
, e in doppia precisione la seconda riga diventa 
perché 
si arrotonda a 
: si ottiene 
e poi 
, completamente sbagliato. Scambiando le righe prima di eliminare il risultato è corretto.
Il pivoting parziale sceglie a ogni passo, come pivot, l’elemento di modulo massimo della colonna, e scambia le righe. Così i moltiplicatori hanno modulo al massimo 1 e il metodo è stabile nella pratica. La fattorizzazione diventa 
, con 
matrice di permutazione. È il metodo standard di tutte le librerie numeriche (LAPACK, NumPy, MATLAB).
Matrici simmetriche definite positive: Cholesky
Se 
è simmetrica e definita positiva (come le matrici di rigidezza delle strutture, le matrici di covarianza, le matrici dei nodi dei circuiti fatti solo di resistori), vale
![\[ A = L\,L^T \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-f3ddffd25bf38f80162273f7cc5ad68b_l3.png "Rendered by QuickLaTeX.com")
con triangolare inferiore. La fattorizzazione di Cholesky costa la metà della LU (
), non richiede pivoting ed è sempre stabile. Se durante il calcolo compare la radice di un numero negativo, la matrice non era definita positiva: è anche il modo più rapido di verificarlo.
Matrici sparse
Nelle applicazioni le matrici grandi sono quasi tutte fatte di zeri: nella matrice dei nodi di un circuito, ogni riga ha elementi non nulli solo per i nodi collegati direttamente. Una matrice 
piena occuperebbe 8 terabyte in doppia precisione; sparsa, con 10 elementi per riga, 80 megabyte. I metodi diretti per matrici sparse riordinano righe e colonne per limitare il riempimento (fill-in), cioè gli zeri che diventano non nulli durante l’eliminazione. È quello che fanno SPICE e i programmi a elementi finiti.
I metodi iterativi
Per matrici enormi e sparse si preferisce spesso partire da una stima 
e migliorarla, senza mai fattorizzare.
Jacobi. Dall’equazione 
si ricava 
usando i valori dell’iterazione precedente per le altre incognite:
![\[ x_i^{(k+1)} = \frac{1}{a_{ii}}\Big(b_i - \sum_{j \ne i}a_{ij}\,x_j^{(k)}\Big) \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-91cf97e1c756c2230a8e90e218f6596d_l3.png "Rendered by QuickLaTeX.com")
Gauss-Seidel. Come Jacobi, ma usa subito i valori già aggiornati nella stessa iterazione. Di solito converge circa due volte più in fretta e usa meno memoria.
Esempio. 
, 
(soluzione 
, 
). Gauss-Seidel da 
: 
, 
; poi 
, 
; poi 
, 
. In tre passi, tre cifre corrette.
Quando convergono
Scrivendo l’iterazione come 
, l’errore viene moltiplicato per 
a ogni passo, quindi il metodo converge se e solo se il raggio spettrale di (il massimo modulo dei suoi autovalori, lezione di Geometria e Algebra) è minore di 1. Una condizione sufficiente facile da verificare: la matrice è a diagonale dominante, 
per ogni riga. Le matrici dei nodi dei circuiti resistivi lo sono quasi sempre.
Per i problemi più grandi si usano metodi iterativi più sofisticati, come il gradiente coniugato per le matrici simmetriche definite positive, che in aritmetica esatta arriva alla soluzione in al massimo 
passi e nella pratica ne richiede molti meno, soprattutto con un buon precondizionatore.
Il residuo non basta
Dopo aver calcolato 
, il residuo 
è facile da calcolare, ma un residuo piccolo non garantisce una soluzione accurata. L’errore relativo è limitato da
![\[ \frac{\|\mathbf x - \hat{\mathbf x}\|}{\|\mathbf x\|} \le K(A)\,\frac{\|\mathbf r\|}{\|\mathbf b\|} \]](https://www.electronics-engineering.com/wp-content/ql-cache/quicklatex.com-8f798c9ead25f00256c9cb0c8a01a299_l3.png "Rendered by QuickLaTeX.com")
con 
il numero di condizionamento della lezione sugli errori. Se è grande, anche un residuo minuscolo può nascondere un errore grande.
