Matrici: Metodi iterativi

I metodi iterativi sono importanti perché spesso risultano più veloci e precisi dei metodi diretti, specialmente per matrici di grandi dimensioni. Essi possono anche essere usati per ``raffinare'' una soluzione ottenuta con un metodo diretto. Inoltre si prestano bene a sfruttare eventuali proprietà di sparsità delle matrici.

Vogliamo trovare ${\bf x}$ tale che $A{\bf x}={\bf b}$ con un metodo iterativo. A tale scopo, scriviamo l'i-esima riga di questo sistema lineare:

$$\sum_{j=1}^n a_{ij}x_j = b_i$$

ed isoliamo a sinistra il termine x_i

$$a_{ii} x_i = b_i - \sum_{j\neq i}a_{ij}x_j .$$ (1)

In un metodo iterativo si parte da una stima `attuale' ${\bf x}$ che in genere non soddisfa quest'ultima equazione. Tuttavia possiamo pensare di trovare una stima migliore ${\bf x}'$ imponendo che l'uguaglianza (1) sia verificata con ${\bf x}'$ a sinistra e ${\bf x}$ a destra; cioè poniamo

$$x_i' = {1\over a_{ii}} \left( b_i -\sum_{j\neq i}a_{ij}x_j \right) .$$

Questo modo di calcolare la nuova approssimazione ${\bf x}'$ è il noto metodo di Jacobi. Per programmarlo basta partire da una stima arbitraria ${\bf x}$ ed aggiornarla via via in base alla relazione data. A tale scopo serve un vettore che memorizzi il ``vecchio'' valore di ${\bf x}$ mentre questo viene aggiornato. Ogni iterazione è dunque del tipo

    for (i=0; i<n; i++) {
      sum = b[i];
      for (j=0; j<n; j++)
          if (j!=i) sum -= a[i][j]*oldx[j];
      x[i] = sum / a[i][i];
    }
    for (i=0; i<n; i++) oldx[i] = x[i];

e questo va inserito in un ciclo while che controlli la convergenza del metodo.

Nel ciclo precedente, è chiaro che quando si calcola x[i] sono già stati calcolati i nuovi valori di x[0],...,x[i-1] ed è possibile usarli nell'aggiornamento di x[i], piuttosto che usare i ``vecchi'' valori in oldx[]. Questo porta al metodo di Gauss-Seidel. Tra l'altro, per programmarlo, si risparmia il vettore oldx, che non serve più.

    for (i=0; i<n; i++) {
      sum = b[i];
      for (j=0; j<n; j++)
          if (j!=i) sum -= a[i][j]*x[j];
      x[i] = sum / a[i][i];
    }

Infine, questo metodo può essere facilmente modificato nel metodo detto ``di rilassamento'' inserendo un peso $\omega$ nell'aggiornamento di x[i]:

      x[i] = (1-w)*x[i] + w*sum/a[i][i];

Scegliendo opportunamente il parametro $\omega$, questo metodo risulta molto più veloce di entrambi i metodi precedenti.

In qualunque metodo iterativo bisogna inserire metodi automatici che controllino la convergenza. Alcune tecniche sono:

• un contatore che fa bloccare il metodo (con fallimento) se si è superato un certo numero prefissato di iterazioni;

• una routine che stima l'errore commesso, ad esempio calcolando $\Vert{\bf x}-{\bf x}'\Vert _\infty$ ad ogni passo, e blocca il metodo (con successo) se questa quantità è inferiore ad un fissato $\varepsilon$;

• tecniche più sofisticate che cercano di stimare il raggio spettrale delle matrici di iterazione per avere dei bound sull'errore.

Si rimanda al libro di testo per lo studio teorico della stima dell'errore.

Una semplice implemetazione di questi metodi iterativi è mostrata nel sorgente iterativi.c. Si noti che per l'implementazione non è richiesta assolutamente la costruzione esplicita delle matrici di iterazione J e G. Piuttosto, bisogna effettuare in modo implicito la moltiplicazione di esse per un vettore. La costruzione esplicita può servire solo per studiarne il raggio spettrale ai fini dello studio teorico del metodo (stimare la convergenza, etc.). Vedere ad esempio il sorgente raggiospettrale.c, che usa il metodo delle potenze per il calcolo dell'autovalore di valore massimo.