Risoluzione numerica dei sistemi lineari: Il metodo del Gradiente Coniugato

 5
direttamente, ma solo tramite la sua applicazione ad un vettore, oppure 
quando non è possibile conservarne in memoria la fattorizzazione. 
La tesi verte sullo studio del metodo del gradiente coniugato, il quale viene 
utilizzato come metodo iterativo, ma, come vedremo, è da considerarsi un 
metodo diretto. 
L’obiettivo è quello di costruire un metodo per la risoluzione di sistemi 
lineari con matrice dei coefficienti simmetrica e definita positiva, chiarire i 
fondamenti matematici alla base del metodo, analizzarne le proprietà di 
convergenza, stabilità, accuratezza e complessità algoritmica, ed illustrare, 
anche attraverso esempi, i vantaggi e i punti deboli. 
Il punto di partenza è l’equivalenza, nel caso di matrice simmetrica definita 
positiva, tra la risoluzione del sistema lineare e la minimizzazione di una 
particolare forma quadratica; allora, il problema è ricondotto a determinare 
il punto di minimo della forma quadratica, partendo da un punto iniziale e 
scegliendo poi opportune direzioni di discesa, lungo le quali muoversi per 
avvicinarsi, il più rapidamente possibile, alla soluzione. La direzione 
ottimale non è ovviamente nota a priori; assumendo particolari direzioni di 
discesa, si costruisce prima il metodo della massima discesa (steepest 
descent), poi il metodo delle direzioni coniugate e da questo, operando una 
scelta più accurata delle direzioni di discesa, si ottiene il metodo del 
gradiente coniugato.  
 6
Si dimostra che, in assenza di errori di round-off, il metodo del gradiente 
coniugato può essere interpretato come un metodo diretto, in quanto 
termina in un numero finito di passi. Tuttavia per matrici di grandi 
dimensioni, esso è generalmente usato come un metodo iterativo perché 
fornisce un’approssimazione accurata della soluzione in un numero di 
iterazioni molto inferiore alla dimensione del sistema. 
Il lavoro comprende, inoltre, l’implementazione in Fortran 90 degli 
algoritmi del metodo di Jacobi e del metodo del gradiente coniugato, con i 
relativi esempi e il confronto tra i risultati. 
La tesi è organizzata come segue:   
 ι  Capitolo 1. Breve introduzione alla risoluzione dei sistemi lineari, ai 
metodi iterativi e al gradiente coniugato. 
 ι  Capitolo 2. Forma quadratica associata ad un sistema lineare. 
 ι  Capitolo 3. Derivazione del metodo della massima discesa. 
 ι  Capitolo 4. Definizioni di autovettori ed autovalori con proprietà e 
derivazione del metodo di Jacobi, con lo studio della convergenza. 
 ι  Capitolo 5. Analisi della convergenza del metodo della massima 
discesa, dagli esempi al caso generale. 
 ι  Capitolo 6. Il metodo delle direzioni coniugate: derivazione, 
convergenza, processo di ortogonalizzazione di Gram-Schmidt. 
 
 7
 ι  Capitolo 7. Derivazione del metodo del gradiente coniugato, con 
analisi della convergenza e della complessità, tecnica del 
precondizionamento ed estensione del metodo al caso di matrice 
qualunque. 
 ι  Capitolo 8. Introduzione al metodo del gradiente coniugato non 
lineare e precondizionamento. 
 ι  Appendice A. Breve nota storica relativa al metodo del gradiente 
coniugato. 
 ι  Appendice B. Vengono presentati gli algoritmi discussi nel lavoro di 
tesi. 
 ι  Appendice C. Esperienze numeriche: implementazioni degli algoritmi 
di Jacobi e del gradiente coniugato, con relativi esempi e confronto 
tra i risultati 
 
 8
CAPITOLO 1 
 
Introduzione alla risoluzione numerica dei sistemi 
di equazioni lineari 
 
 
1.1     L’algebra lineare numerica e i sistemi lineari 
L’obiettivo della Matematica Computazionale è quello di fornire strumenti 
software per la risoluzione effettiva di problemi concreti. Nel processo di 
risoluzione di un problema gioca un ruolo fondamentale la fase di 
risoluzione numerica del modello matematico che descrive il problema in 
esame, la fase che va dalla formulazione del modello numerico fino al 
sviluppo del software. I più semplici modelli matematici sono di tipo 
lineare: essi sono i più utilizzati e sono di fondamentale importanza in 
quanto l’analisi e lo studio di problemi concreti, anche estremamente 
complessi, sono in primo luogo effettuati su un modello lineare 
semplificato. 
L’Algebra Lineare Numerica ha come obiettivo lo studio, lo sviluppo e 
l’analisi di metodi computazionali e delle loro tecniche implementative per 
la risoluzione di problemi concreti della scienza e della tecnica descritti 
mediante modelli matematici lineari. Tali problemi derivano direttamente 
 9
dall’applicazione di leggi che regolano il fenomeno in esame (come accade 
in molti campi della scienza, quali fisica, biologia, economia) oppure 
compaiono nel processo di risoluzione di problemi di altri settori 
dell’Analisi Numerica, quali interpolazione, approssimazione, 
ottimizzazione, risoluzione di problemi differenziali. 
E’ quindi evidente il ruolo predominante dell’Algebra Lineare Numerica 
nell’ambito della Matematica Computazionale. 
I problemi fondamentali dell’ Algebra Lineare Numerica sono i seguenti : 
 ξ  risoluzione di sistemi di equazioni lineari 
 ξ  risoluzione di problemi di autovalori ed autovettori 
 ξ  risoluzione di problemi di minimi quadrati 
Il sistema lineare è uno strumento di base della modellistica matematica. 
Molti problemi scientifici, infatti, sono descritti mediante un modello 
matematico basato su un sistema di equazioni lineari 
Il processo di risoluzione di tali tipi di problemi consiste, in generale, di tre 
passi: 
1. formulare il modello matematico del problema                                       
(esprimere in termini matematici le relazioni esistenti tra le quantità, 
note ed incognite, che intervengono nel problema). 
2. calcolare le quantità note del problema. 
3. risolvere un sistema di equazioni lineari per ottenere la soluzione del 
problema. 
 10
L’esecuzione dei primi due passi dipende dalla conoscenza del problema, 
mentre l’esecuzione del terzo passo consiste nella scelta e nell’utilizzo di un 
opportuno algoritmo e del relativo software. 
 
 11
1.2 Risoluzione numerica dei sistemi lineari : Metodi Iterativi  
La risoluzione dei sistemi di equazioni lineari è uno dei problemi 
dell’Analisi Numerica che si presenta più spesso nelle applicazioni. La 
scelta dell’algoritmo risolutivo più efficace richiede un’analisi delle 
caratteristiche del sistema, ad esempio numero e tipo delle equazioni. In 
molte applicazioni, le caratteristiche del modello numerico utilizzato per 
descrivere il problema da risolvere e la necessità di ottenere risultati 
accurati, possono condurre a sistemi lineari di grandi dimensioni, cioè con 
un grande numero di equazioni e di incognite, e sparsi, cioè i coefficienti di 
molte incognite sono nulli e quindi ogni equazione coinvolge solo un 
piccolo numero di incognite. 
I metodi risolutivi per i sistemi di equazioni lineari si distinguono in due 
tipi: metodi diretti e metodi iterativi. Nei metodi diretti si giunge alla 
soluzione esatta (a meno degli errori di arrotondamento) con un numero 
finito di operazioni sui dati; nei metodi iterativi la soluzione viene invece 
approssimata dai termini di una successione di cui la soluzione cercata è il 
limite. La convenienza dell’uno o dell’altro tipo di metodo dipende da 
particolari proprietà del sistema che si vuole risolvere. 
Se il sistema lineare è di grandi dimensioni, può risultare poco efficiente o 
addirittura impraticabile la sua risoluzione mediante un metodo diretto, 
come ad esempio il metodo di Gauss che consiste nel trasformare il sistema 
di partenza in un sistema triangolare equivalente, mediante il metodo di 
 12
eliminazione di Gauss, e nel risolvere successivamente il sistema 
triangolare mediante il metodo di back-substitution. Infatti, il metodo di 
Gauss richiede un numero di operazioni floating-point tale che per sistemi 
di ordine n >10
4
, il tempo totale di CPU relativo alla risoluzione del sistema 
su un calcolatore può anche essere dell’ordine di qualche giorno. 
Inoltre, se il sistema è sparso, quindi la matrice dei coefficienti è sparsa, 
può risultare svantaggioso risolvere il sistema con un metodo diretto. Se a 
un tale sistema si applica un metodo diretto, le matrici dei sistemi intermedi 
o di arrivo possono diventare matrici dense, cioè con un elevato numero di 
elementi non nulli. Ad esempio, il metodo di Gauss durante il processo di 
trasformazione del sistema dato nel sistema triangolare equivalente, 
modifica i coefficienti del sistema ed a causa di tali modifiche può accadere 
che coefficienti uguali a zero diventino non nulli (fenomeno detto fill-in). 
 Sorgono così seri problemi di costo computazionale e di ingombro di 
memoria ed è allora importante utilizzare metodi che sfruttino, dal punto di 
vista della complessità computazionale, il fatto che la maggior parte dei 
coefficienti sono nulli.  
Una valida alternativa ai metodi diretti, nella risoluzione di sistemi lineari 
con matrice dei coefficienti di ordine elevato e sparsa, è costituita dai 
metodi iterativi. Essi hanno il vantaggio di preservare la struttura della 
matrice, sfruttandone quindi meglio l’eventuale sparsità, e di essere, in generale, 
di più semplice implementazione. 
 13
I metodi iterativi per la risoluzione di un sistema lineare  
bAx   
generano a partire da un vettore iniziale 
)0(
x , una successione di vettori 
{
)(k
x }che, sotto opportune ipotesi, converge alla soluzione del problema. 
Nei metodi diretti, che valutano la soluzione in un numero finito di passi, 
gli errori presenti nei risultati nascono esclusivamente dall’aritmetica finita 
e/o dalla presenza di errori sui dati. In quelli iterativi, invece, agli errori 
sperimentali e di arrotondamento si aggiungono gli errori di troncamento, 
derivanti dal fatto che il limite cercato deve essere necessariamente 
approssimato troncando la successione per un indice sufficientemente 
grande. Si aggiunge, quindi, il problema della determinazione di un criterio 
di arresto. 
I metodi iterativi risultano, dunque, convenienti per matrici di grandi 
dimensioni, specialmente se sparse. Essi non richiedono una modifica della 
matrice del sistema e nemmeno la sua effettiva memorizzazione; è solo 
necessario poter accedere in qualche modo ai suoi elementi. Inoltre, i 
metodi iterativi possono essere utilizzati nei casi in cui si voglia raffinare 
una soluzione approssimata ottenuta con altri algoritmi, o mediante 
informazioni a priori sul problema. Infine, è possibile ridurre di molto il 
tempo di elaborazione, eseguendo un minor numero di iterazioni in quei 
casi in cui non sia richiesta un’elevata accuratezza. 
 14
1.3     Il metodo del gradiente coniugato 
Il metodo del gradiente coniugato è il più importante metodo iterativo per 
la risoluzione di sistemi lineari sparsi e di grandi dimensioni. Esso è 
applicato a sistemi della forma 
bAx   
dove x  è un vettore incognito, b  è un vettore noto ed A  è una matrice nota, 
quadrata, simmetrica e definita positiva. Questi sistemi si presentano in 
molte importanti applicazioni, come i metodi delle differenze finite e degli 
elementi finiti per la risoluzione approssimata di problemi ai limiti per 
equazioni differenziali a derivate parziali, l’analisi strutturale e l’analisi dei 
circuiti. 
Si può dimostrare che, in aritmetica infinita, per un sistema lineare di ordine 
n ,  il metodo del gradiente coniugato converge alla soluzione esatta in al 
più n  iterazioni ed è quindi da considerare un metodo diretto. Viene però di 
fatto utilizzato come un metodo iterativo perché, specie se opportunamente 
precondizionato, fornisce un’approssimazione sufficientemente accurata 
della soluzione in un numero di iterazioni molto inferiore ad n  (ad ogni 
iterazione fornisce un’approssimazione migliore della soluzione). Ciò è 
particolarmente utile quando si opera su sistemi lineari di dimensione 
estremamente elevata, per i quali sarebbe inaccettabile effettuare un numero 
di iterazioni pari al numero delle equazioni.  
 15
Va inoltre osservato che, a causa degli errori provocati dall’aritmetica di 
macchina, il metodo non riuscirebbe comunque a produrre la soluzione 
esatta in n  iterazioni. 
 
 16
CAPITOLO 2 
 
La forma quadratica 
 
 
2.1 Notazioni e definizioni 
Consideriamo il sistema lineare 
                                                       bAx                                                     (1) 
dove A  è una matrice nn υ, x  e b  sono vettori di dimensione n                  
( cioè matrici 1 υn  ). 
Scrivendo il sistema per esteso, l’equazione (1) diventa 
                                            
 ÷
 ÷
 ÷
 ≠
 •
 ♦
♦
 ♦
 ♥
 ♣
nnn
n
aa
aa
...
::
...
1
111
 
 ÷
 ÷
 ÷
 ≠
 •
 ♦
♦
 ♦
 ♥
 ♣
n
x
x
:
1
  =  
 ÷
 ÷
 ÷
 ≠
 •
 ♦
♦
 ♦
 ♥
 ♣
n
b
b
:
1
 
Dati due vettori x  e y , si denota con yx
T
  il prodotto scalare (o interno) di 
x  e y ,  cioè la somma scalare: 
 ƒ
  
  
n
i
ii
T
yxyx
1
 
Notiamo che  xyyx
TT
    e se x  e y  sono ortogonali, allora 0  yx
T
. 
 
 17
Data una matrice quadrata A  di ordine n ,  
1. si denota con 
T
A  la trasposta di A , cioè la matrice 
T
AC   tale che : 
jiij
ac   . 
2. A   è simmetrica se  
T
AA  . 
3. A   è definita positiva se, per ogni vettore non nullo x , si ha :  
                                                   Axx
T
> 0                                                 (2) 
 
 18
2.2 La forma quadratica 
Una forma quadratica è una funzione scalare f , di secondo grado, di un 
vettore x , espressa nella forma : 
 
                                         cxbAxxxf
TT
      
2
1
)(                                    (3) 
 
dove A  è una matrice, x  e b  sono vettori, e c  è una costante.  
 
Dimostreremo che se A  è simmetrica e definita positiva, la risoluzione del 
sistema lineare bAx   è equivalente a trovare il punto di minimo della 
forma quadratica )(xf  (detta anche energia del sistema). 
 
Esempio : 
                                  0,
8
2
,
62
23
  
 ÷
 ÷
 ≠
 •
 ♦
 ♦
 ♥
 ♣
  
  
 ÷
 ÷
 ≠
 •
 ♦
 ♦
 ♥
 ♣
  cbA                                      (4)