Eliminazione di Gauss

Il metodo di Gauss è uno dei più semplici e più usati (nella variante col pivot). Esso costruisce la fattorizzazione LUmoltiplicando successivamente A a sinistra per delle matrici elementari di Gauss. Ciascuna moltiplicazione annulla gli elementi che si trovano sotto il pivot, attraverso una combinazione lineare delle righe di A. Questo processo corrisponde ad espandere Lnella forma $L=M_1^{-1}\cdots M_n^{-1}$, e moltiplicare ad ogni passo per un certo M_i. Vedi pagina 157 del libro di testo.

Al generico passo k vengono annullati gli elementi sotto-diagonali della k-esima colonna: questo corrisponde ad eliminare la variabile x_k dalle corrispondenti equazioni, da cui il nome del metodo.

Prima di compiere il k-esimo passo la matrice A è stata trasformata nella forma

$$\left[ \matrix{U_{k-1} & A' \cr 0 &A''} \right]$$

con U_k-1 triangolare superiore. Il k-esimo passo consiste nel moltiplicare a sinistra per una matrice elementare M_k che differisce da I solo per gli elementi sotto-diagonali della k-esima colonna, che sono

$$-{a_{k+1,k}\over a_{kk}}, -{a_{k+2,k}\over a_{kk}}, \ldots, -{a_{nk}\over a_{kk}}.$$

La moltiplicazione dunque consiste nel sostituire la i-esima riga di A con una combinazione di essa con la k-esima, tale da permettere di annullare l'elemento a_ik. Si può verificare che le prime k righe di A rimangono invariate (perché le prime k righe di M_k sono uguali all'identità). Indicando con A_(i) la i-esima riga di A, l'operazione è:

$$A_{(i)} \leftarrow A_{(i)} - {a_{ik}\over a_{kk}} A_{(k)} \hspace{2cm}\mbox{ per $i>k$}.$$

Osserviamo che, per ogni riga i, abbiamo

$$a_{ik} \leftarrow a_{ik} - {a_{ik}\over a_{kk}} a_{kk} = 0 \hspace{2cm}\mbox{ per $i>k$},$$

cioè tutti gli elementi sotto a_kk sono stati annullati.

A causa della struttura di A, anche le prime k-1 colonne di A rimangono invariate, quindi l'aggiornamento va effettivamente eseguito solo sulle colonne di indice $\geq k$. Il frammento di programma che effettua il passo k sulla matrice Aè dunque

    /* Passo k dell'eliminazione di Gauss */
    for (i=k+1; i<n; i++)
        for (j=k; j<n; j++)
            a[i][j] = a[i][j] - a[i][k]*a[k][j]/a[k][k];

Per evitare inutili accessi al vettore a è meglio porre in variabili locali alcune quantità che rimangono costanti all'interno di ciasun ciclo for, in particolare il pivot a_kked il valore m=a_ik/a_kk. Otteniamo:

    /* Passo k dell'eliminazione di Gauss */
    pivot = a[k][k];
    for (i=k+1; i<n; i++) {
        m = -a[i][k] / pivot;
        for (j=k; j<n; j++)
            a[i][j] = a[i][j] + m*a[k][j];
    }

L'intero processo di eliminazione consiste dunque nelle precedenti righe, iterate dentro un ciclo for (k=0; k<n; k++). All'uscita di questo ciclo A sarà stata trasformata nella matrice triangolare superiore U.

I valori assunti da m nell'ultimo frammento sono gli elementi della matrice L (a meno del segno). Se li si vuole memorizzare, si può scrivere, dentro il ciclo su i:

        l[i][k] = a[i][k] / pivot;

Visto però che A alla fine del calcolo sarà triangolare superiore (diverrà uguale ad U), mentre L è triangolare inferiore, gli elementi di L possono essere memorizzati nelle posizioni [i][k] della stessa matrice A. Così, alla fine, nelle locazioni dove prima era Asi troverà la matrice L (a cui manca però la diagonale di uni) accostata alla matrice U. In questo modo si ottiene la fattorizzazione LU di A. Per l'implementazione vedere il sorgente lu.c.