Per proseguire il discorso iniziato il mese scorso, vedremo ora i metodi iterativi che, come già si è detto, sono utilizzati quando la matrice del sistema Ax=b è sparsa, di grandi dimensioni o priva di una struttura particolare. La soluzione del problema è ottenuta come limite di una successione di soluzioni di sistemi più semplici.
Ciò che sta alla base dei metodi iterativi è lo splitting della matrice A: in sostanza si considera A=P-N, dove P è una matrice non singolare. Da ciò segue che se Ax=b è il sistema dato,
sarà la forma che ci permetterà di esprimere la successione delle soluzioni in forma ricorsiva; infatti, dato il vettore x0, detto stima iniziale, si ricava xk per k”1 risolvendo il seguente sistema:
Ottenuta in questo modo la successione {xk}, si ha che
dove x è la soluzione del sistema Ax=b. In particolare, dalla definizione di limite di una successione si ha che la (3) è verificata quando xk-x<®. Ponendo, infine, ek=xk-x si dice che il metodo è convergente quando ek<®.
Vediamo ora una condizione sufficiente per la convergenza.
Riscrivendo la (1) nella seguente forma
e sottraendola dalla (2) si ottiene
da cui, ponendo B=P-1N, dove B è la matrice di iterazione del metodo, risulta
Iterando la (6) si ottiene la relazione
dalla quale, sapendo che e0 è il vettore costante che rappresenta l’errore commesso nel dare la stima iniziale della soluzione
(e0=x0-x), si ricava che:
lim(k”‘)ek+1=
=lim(k”‘)Bk+1e0=
=e0lim(k”‘)Bk+1
dunque lim(k”‘)ek+1=0 se e soltanto se
lim(k”‘)Bk+1=0.
Questo teorema, però, ci permette di valutare la convergenza del metodo studiando una successione a valori nello spazio n-dimensionale; tuttavia, essendo più semplice lavorare nello spazio reale unidimensionale, può essere utile sapere che:
a. lim(k”‘)Bk=0
b. «(B)<1
c. lim(k”‘)Bk+1=0 per qualche
sono equivalenti. Sfruttando dunque le proprietà (b) e (c) si può ottenere più agevolmente una stima della convergenza di B; infatti, non occorre più valutare il limite di una successione a valori nello spazio reale n-dimensionale ma, nel primo caso, si valuta una disuguaglianza nello spazio reale e nel secondo il limite di una successione dello spazio unidimensionale.
Altra caratteristica della succesione {ek} è che la sua convergenza è tanto più rapida quanto più piccolo è «(B).
Ritornando allo splitting della matrice A occorre specificare che a seconda di come P ed N vengono scelte, si ottengono differenti metodi iterativi, l’applicazione dei quali al sistema Ax=b può portare talora a differenti risultati a causa degli eventuali errori prodotti dal metodo sulla matrice specifica. Supponendo di suddividere gli elementi di A come indicato in figura, vediamo quali possono essere le scelte di P ed N ed i metodi corrispondenti:
1. P=D N=-(E+F)
– Metodo di Jacobi
– La componente i-esima dell’iterata
(k+1)-esima è data dalla seguente
relazione
ottenibile esplicitando
xk+1=P-1Nxk+P-1b.
– La matrice di iterazione del metodo è
BJ=-D-1(E+F)
– E’ possibile dimostrare che:
condizione sufficiente, ma non necessaria, per la convergenza del metodo è che A sia a dominanza diagonale stretta.
– Questo metodo è anche detto metodo degli spostamenti simultanei in quanto per il calcolo della i-esima componente del vettore xk+1 occorrono tutte le componenti del vettore precedente tranne l’i-esima.
Questo implica che solo dopo aver calcolato tutte le componenti di xk+1 è possibile liberare l’area di memoria contenente xk e l’operazione viene compiuta contemporanemente su tutte le sue componenti.
2. P=D+E N=-F
– Metodo di Gauss-Sidel
– La componente i-esima dell’iterata (k+1)-esima è data dalla seguente relazione
– La matrice di iterazione del metodo è BGS=-(D+E)-1F.
– E’ possibile dimostrare che:
a. se A è a dominanza diagonale stretta il metodo converge
b. se A è simmetrica e definita positiva il metodo converge (questo, in genere, non si verifica per il metodo di Jacobi)
c. se A è tridiagonale si ha che «(BGS)=«2(BGS) quindi in generale il metodo di Gauss-Sidel converge più velocemente del metodo di Jacobi
– Questo metodo è anche detto metodo degli spostamenti successivi in quanto, per il calcolo di xk+1i, servono le prime i-1 componenti di xk
(si noti che non si utilizzano xk1,…,xkn); questo permette di operare utilizzando un solo vettore nel quale si sostituisce successivamente l’i-esima componente di xk con l’i-esima componente di xk+1.
