Analisi numerica/Sistemi lineari

Nel seguito, salvo indicazioni contrarie, sia $A$ una matrice quadrata di ordine $n$, ${\underline {b}}$ e ${\underline {x}}$ due vettori di $\mathbb {R} ^{n}$. Chiameremo sistema lineare di n incognite in n equazioni (o sistema lineare di ordine n) una relazione del tipo

$A{\underline {x}}={\underline {b}}.$

Tale sistema, come è noto, ammette un'unica soluzione se e solo se $det(A)$$\neq$$0$; in caso contrario la soluzione esiste (non unica) solo se ${\underline {b}}$$\in$$range(A)$. Nel seguito ci riferiremo principalmente a sistemi in cui la matrice abbia determinante non nullo.

Metodi diretti

In questa sezione analizzeremo i metodi diretti, ossia quei metodi che non costruiscono una successione di approssimanti che converga verso la soluzione del sistema, ma cercano piuttosto di effettuare una serie di passaggi algebrici per giungere direttamente alla soluzione. In genere si tratta di fattorizzazioni della matrice, ossia delle scomposizioni che tendono a ricondurre il sistema a due o più sistemi da risolvere in successione, la cui complessità sia però molto inferiore a quella del sistema di partenza.

Sistemi triangolari

I sistemi triangolari sono quelli in cui la matrice abbia la struttura triangolare (inferiore o superiore). Essi sono particolarmente facili da risolvere, anche a mano, se si utilizza l'algoritmo delle sostituzioni successive. Vediamo ad esempio come si risolve un sistema triangolare superiore (il caso inferiore è completamente speculare e il lettore può provare a impostarne la risoluzione come esercizio). Sia dato dunque un sistema della forma

${\begin{pmatrix}a_{1,1}&a_{1,2}&\dots &a_{1,n}\\0&a_{2,2}&\dots &a_{2,n}\\\vdots &\vdots &\ddots &\vdots \\0&0&\dots &a_{m,n}\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\\\vdots \\x_{n}\end{pmatrix}}={\begin{pmatrix}b_{1}\\b_{2}\\\vdots \\b_{n}\end{pmatrix}}$

Appare evidente che l'ultima equazione si risolve rispetto a $x_{n}$ con una sola operazione, dividendo il termine noto per il coefficiente $a_{nn}$. Ma fatto ciò, nella penultima equazione il termine ${a_{(n-1)n}}x_{n}$ può essere portato a secondo membro, in quanto completamente noto, e si può risolvere come prima rispetto a $x_{n-1}$, e così via. Il costo computazionale di tale risoluzione è abbastanza basso: all'i-esima equazione (partendo dal basso) occorre fare $i-1$ moltiplicazioni, $i-1$ somme e una divisione, per un totale di $2i-1$ operazioni. Sommando per i che va da 1 a $n$, otteniamo

$\#flops=2{\frac {n(n+1)}{2}}-n=n^{2}.$

Il metodo di eliminazione gaussiana

Data la facilità con cui si risolve un sistema triangolare, cercare di ricondurre il sistema ad una serie di sistemi triangolari sembrerebbe una buona idea. Ovviamente tale procedimento non deve essere troppo costoso, altrimenti la strategia perde di utilità. Un metodo particolarmente adatto a questo scopo è il cosiddetto metodo di eliminazione gaussiana (in breve MEG), che si impara a fare a mano, per piccoli sistemi, già alle scuole superiori. La strategia è molto semplice e si basa su due semplici proprietà delle equazioni:

• moltiplicando entrambi i membri di un'equazione per un numero diverso da zero, la soluzione non cambia.
• sommando a entrambi i membri di un'equazione una stessa quantità, la soluzione non cambia.

Sfruttando queste due semplici proprietà, si cerca di ridurre la complessità delle equazioni in modo successivo: si comincia dalla prima equazione, esprimendo $x_{1}$ in funzione delle altre $n-1$ incognite, dopodiché si sostituisce la sua espressione nelle altre equazioni. In questo modo si ottiene un sotto-sistema di ordine $n-1$, cui si può applicare lo stesso procedimento, fino ad ottenere un sistema banale di ordine 1, che si risolve in un passaggio. Fatto questo, sostituendo a ritroso, si ottiene il valore di tutte le altre incognite. Vediamo di fare un esempio pratico, per capire come opera a livello algebrico il MEG. Supponiamo di voler risolvere il sistema

${\begin{pmatrix}3&1&2\\1&1&1\\1&0&2\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\\x_{3}\end{pmatrix}}={\begin{pmatrix}1\\0\\2\end{pmatrix}}$

La prima cosa che vogliamo fare è di eliminare $x_{1}$ dalla seconda e terza equazione. Per fare questo, dobbiamo combinarle linearmente con la prima in modo da annullare il coefficiente della prima colonna. Se sommiamo la seconda alla prima moltiplicata per $-a_{2,1}/a_{1,1}$ e sommiamo la terza alla prima moltiplicata questa volta per $-a_{3,1}/a_{1,1}$, otteniamo un nuovo sistema della forma

${\begin{pmatrix}3&1&2\\0&2&1\\0&-1&4\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\\x_{3}\end{pmatrix}}={\begin{pmatrix}1\\-1\\5\end{pmatrix}}$

Come si può osservare, le operazioni che vengono fatte sulle righe della matrice vanno fatte anche sul termine noto. A questo punto la seconda e terza equazione sono indipendenti da $x_{1}$. Ripetiamo il procedimento per eliminare $x_{2}$ dalla terza equazione, sommandola alla seconda moltiplicata per $-a_{3,2}/a_{2,2}$, ottenendo

${\begin{pmatrix}3&1&2\\0&2&1\\0&0&9/2\end{pmatrix}}{\begin{pmatrix}x_{1}\\x_{2}\\x_{3}\end{pmatrix}}={\begin{pmatrix}1\\-1\\9/2\end{pmatrix}}$

Ci siamo dunque ricondotti ad un sistema triangolare superiore, la cui risoluzione richiede $n^{2}$ operazioni con l'algoritmo delle sostituzioni successive e il problema può dunque considerarsi risolto.

Fattorizzazioni

Introduciamo ora i metodi di risoluzione diretta, i quali operano sostanzialmente fattorizzando la matrice, ossia esprimendola come prodotto di matrici dotate di strutture particolari, in modo da ridurre il sistema alla risoluzione di più sistemi, ma di complessità decisamente inferiore.

Fattorizzazione LU

Anche il MEG ha un'interpretazione fattorizzata. Più precisamente ci si riferisce al MEG come metodo di fattorizzazione LU, dove le due lettere stanno per Lower e Upper, in quanto il MEG coincide con un metodo di fattorizzazione in due matrici, la prima triangolare inferiore e la seconda triangolare superiore. Per capire come sono fatte le matrici del MEG, ripensiamo ai passi dell'algoritmo.

All'i-esimo passo, dobbiamo eliminare la variabile $x_{i}$ dalle equazioni successive. Per far questo prendiamo la generica riga j-esima (con $j>i$) e la sommiamo alla prima riga, premoltiplicata per $-a_{j,i}/a_{i,i}$. Quindi stiamo premoltiplicando la matrice per il vettore riga

${\begin{pmatrix}\cdots &0&-a_{j,i}/a_{i,i}&0&\cdots &0&1&0&\cdots \end{pmatrix}}.$

In realtà, volendo eliminare tutte le variabili con indice minore di j dalla riga j-esima in un'unica operazione, dobbiamo premoltiplicare per la matrice

$M_{j}={\begin{pmatrix}1&\cdots &0&0&\cdots &0\\\vdots &\ddots &\vdots &\vdots &&\vdots \\0&&1&0&&0\\0&&-a_{j+1,j}/a_{j,j}&1&&0\\\vdots &&\vdots &\vdots &\ddots &\vdots \\0&\cdots &-a_{n,j}/a_{j,j}&0&\cdots &1\end{pmatrix}}.$

Dunque, definendo la matrice

$M=M_{n}M_{n-1}\cdots M_{2}M_{1},$

il MEG opera la trasformazione $MA=U$, da cui ricaviamo che $L=M^{-1}$. Ora, il calcolo di un'inversa è in genere brutto numericamente, ma data la particolare forma delle $M_{j}$, abbiamo che

$M_{j}^{-1}={\begin{pmatrix}1&\cdots &0&0&\cdots &0\\\vdots &\ddots &\vdots &\vdots &&\vdots \\0&&1&0&&0\\0&&a_{j+1,j}/a_{j,j}&1&&0\\\vdots &&\vdots &\vdots &\ddots &\vdots \\0&\cdots &a_{n,j}/a_{j,j}&0&\cdots &1\end{pmatrix}},$

per cui anche la costruzione della matrice $L$ risulta abbastanza semplice.

Una volta completata la fattorizzazione, bisogna risolvere in sequenza

$L{\underline {y}}={\underline {b}}$

$U{\underline {x}}={\underline {y}}$

per ottenere la soluzione del problema di partenza. Il costo computazionale della fattorizzazione è di $O(2/3n^{3})$ flops, cui poi si aggiungono i $2n^{2}$ flops necessari per la risoluzione dei due sistemi triangolari. La maggior parte del costo computazionale risiede quindi nella fattorizzazione. Questo va tenuto in considerazione, nel caso si debbano risolvere più sistemi lineari aventi la stessa matrice, in modo da effettuare una volta sola la fattorizzazione.

Analizziamo alcune proprietà della fattorizzazione LU. Cominciamo con l'esporre il risultato principale.

Data una matrice A, esiste un'unica fattorizzazione LU se e solo se i minori principali di ordine $1,...n-1$ hanno tutti determinante non nullo.

Notiamo che non è richiesto che la matrice A sia non singolare. Se A è singolare, giunti all'ultimo passo della fattorizzazione (ammettendo che ciò sia possibile) ci si troverà con $a_{nn}=0$, ma dato che non ci sono passaggi che richiedono la divisione per tale coefficiente, la fattorizzazione esiste comunque. Ovviamente il sistema lineare non sarà risolvibile in modo univoco; sarà indeterminato se anche $b_{n}$ risulterà nullo, altrimenti sarà impossibile.

Metodi iterativi

Nell'analisi numerica in generale, si dicono iterativi, tutti quei metodi che costruiscono delle successioni di approssimanti della soluzione del problema reale. Salvo casi particolari, dunque, nessuna delle approssimazioni generate coincide con la soluzione esatta del problema; quello che però il metodo iterativo punta a soddisfare è un criterio di convergenza della soluzione approssimata, in una norma appropriata, alla soluzione esatta.

Nel caso dei sistemi lineari, dato che si lavora in $\mathbb {R} ^{n}$, tutte le norme sono equivalenti. Tuttavia le norme più comunemente utilizzate per ottenere stime di convergenza sono quella euclidea e quella dell'energia; quest'ultima è definita solo per sistemi con matrice simmetrica e definita positiva in quanto è data da

$||{\underline {x}}||_{A}=<{\underline {x}},A{\underline {x}}>=||A^{1/2}{\underline {x}}||_{2}$

dove $<\cdot ,\cdot >$ è l'usuale prodotto scalare in $\mathbb {R} ^{n}$, mentre $||\cdot ||_{2}$ è la norma euclidea.