INTRODUZIONE V
sistemi di equazioni differenziali, ordinarie o a derivate parziali.
Se la discretizzazione Ł abbastanza fine, allora essa pu dar luogo
ad un sistema lineare di dimensione molto grande. D'altra parte la
particolare natura delle equazioni differenziali insieme a
particolari scelte delle discretizzazioni, danno luogo
generalmente a sistemi lineari di forma particolare che
consentono l'uso di certe tecniche iterative che sono presentate e
studiate nel capitolo 1 paragrafi 1.4-1.5-1.6-1.7.
Per il fatto che il sistema lineare Ł uno strumento di base della
modellistica matematica, ne scaturisce la necessit di avere a
disposizione algoritmi efficienti e per i quali siano note le
propriet di stabilit . Alla variet di contesti in cui pu nascere
un sistema lineare, corrisponde un altrettanta variet di tipi di
sistemi da risolvere. Per esemplificare, i sistemi che provengono
dalla risoluzione di equazioni alle derivate parziali sono
sostanzialmente diversi, per dimensione, propriet di simmetria o
numero di elementi nulli nella matrice dei coefficienti, da quelli
che si ottengono nei modelli di tipo stocastico per la simulazione
di code. Anche nell ambito di una stessa categoria di problemi,
come quelli relativi alle derivate parziali, il tipo di sistema pu
variare a seconda del metodo di discretizzazione utilizzato:
differenze finite, metodi spettrali o elementi finiti.
Inoltre anche nell ambito di un particolare metodo, come negli
elementi finiti, dal tipo di elemento utilizzato.
E importante, pertanto avere a disposizione una variet di
algoritmi, dai quali poter scegliere il piø adatto, per stabilit ,
INTRODUZIONE VI
occupazione di memoria, velocit di esecuzione, a risolvere su un
determinato calcolatore un sistema particolare. In effetti la
stabilit dell algoritmo, il suo costo e il tipo di calcolatore
utilizzato rappresentano i criteri piø importanti per una scelta
ottimale.
Lo scopo principale di questa tesi Ł la valutazione delle
performances dei metodi iterativi per la risoluzione dei sistemi
lineari.
Abbiamo preso in esame i metodi iterativi di Jacobi e Gauss-
Seidel, la tecnica del rilassamento ed il metodo del gradiente
coniugato.
I primi due generano una successione di soluzioni approssimanti,
a partire da uno split della matrice del sistema linerare.
La tecnica del rilassamento Ł un evoluzione del metodo di
Gauss-Seidel per accelerare la velocit di convergenza.
Il metodo del gradiente coniugato trasforma il problema di
partenza in un problema di ricerca del minimo di un funzionale.
Nei test abbiamo scelto particolari matrici per stressare le
prestazioni dei metodi e per metterne in risalto le propriet .
Per implementare i test abbiamo scritto un codice Visual Basic
che un front-end accattivante permette in maniera amichevole
la risoluzione di un sistema lineare con i metodi iterativi trattati
nella tesi.
CAPITOLO 1
I metodi iterativi
SOMMARIO: 1.1 I sistemi di equazioni lineari.-1.2 Metodi diretti.-1.3 Metodi iterativi.-
1.4.-Metodo di Jacobi.-1.5 Metodo di Gauss-Seidel.-1.6 Metodo di rilassamento (SOR).-1.7
Metodo del gradiente coniugato.
1.1 I sistemi di equazioni lineari
Il problema della risoluzione di sistemi di equazioni lineari si
presenta in moltissime applicazioni, sia esplicitamente nel
modello matematico associato al fenomeno fisico in esame, sia
come passo intermedio o finale nella risoluzione numerica del
modello in questione, rappresentato per esempio da equazioni
differenziali. Sistemi lineari sono anche coinvolti, nella stessa
costruzione di metodi numerici per il calcolo di autovalori di
matrici, nella risoluzione di sistemi di equazioni non lineari, di
equazioni differenziali, nell approssimazione di dati e di
funzioni. Un sistema del tipo:
A x = b con A ∈ Rn x m, b, x ∈ Rn
ammette soluzione se e solo se il vettore b appartiene allo spazio
lineare generato dalle colonne di A.
CAPITOLO 1. I METODI ITERATIVI 2
Infatti, posto:
A = (a1 a2, .an),con ai ∈ Rn e x =(x1, x2, xn)T xi ∈ Rn ,
la scrittura A x= b pu essere interpretata nella forma
b = (x1 a1+ x2 a2+ .+ xn an).
Inoltre, la soluzione Ł unica se e solo se la matrice A Ł non
singolare ossia det (A)≠0, in questo caso possiamo scrivere
x = A-1 b.
Di norma i metodi numerici per la risoluzione dei sistemi lineari
vengono suddivisi in due classi: metodi diretti e metodi iterativi.
Con i metodi diretti l esatta soluzione viene costruita, in assenza
di errori di arrotondamento, in un numero finito di passi e si tratta
sostanzialmente dei metodi che utilizzano l idea
dell eliminazione di Gauss. Per sistemi A x = b con matrici A
dense, cioŁ con la maggior parte degli elementi non nulli, i
metodi diretti sono di solito piø efficienti.
I metodi iterativi invece sono generalmente utilizzati per la
risoluzione di sistemi con matrici A sparse, cioŁ con molti
elementi nulli e di ordine elevato. Sistemi sparsi sono presenti in
numerose applicazioni. A volte, come infatti succede nella
risoluzione numerica di equazioni alle derivate parziali con
metodi alle differenze finite, gli elementi non nulli costituiscono
delle successioni regolari di pochi numeri distinti. A causa
dell elevato ordine delle matrici coinvolte in problemi di questo
tipo i metodi diretti non sempre sono utilizzabili. Infatti, questi
ultimi, ottengono la soluzione x in numero finito di passi
mediante una successione finita di trasformazioni del problema
CAPITOLO 1. I METODI ITERATIVI 3
iniziale in problemi equivalenti, cioŁ con la stessa soluzione x,
ma con matrici dei coefficienti diverse; anzi , col procedere del
metodo di risoluzione il numero di elementi non nulli presenti in
queste matrici generalmente cresce, e pu ben presto saturare lo
spazio disponibile della memoria centrale del calcolatore. In
questi casi Ł utile, e spesso indispensabile, utilizzare metodi
iterativi, i quali, operando sempre e solo con gli elementi della
matrice iniziale A, generano una successione infinita di vettori
convergente, sotto opportune condizioni, alla soluzione cercata.
PoichØ il processo iterativo lascia inalterata la matrice A, Ł
sufficiente memorizzare gli elementi non nulli di A. Osserviamo
infine che, mentre i metodi diretti in assenza di errori nella
rappresentazione dei numeri forniscono la soluzione esatta del
sistema, indipendentemente da quella che potrebbe essere la
precisione richiesta dall utente, i metodi iterativi, se convergenti,
consentono invece, di arrestare il processo iterativo non appena
la precisione desiderata Ł stata raggiunta. Va tuttavia sottolineato
che un metodo iterativo risulter efficiente solo se consentir di
raggiungere la precisione desiderata con un numero accettabile di
iterazioni ovvero di operazioni aritmetiche.
CAPITOLO 1. I METODI ITERATIVI 4
1.2 Metodi diretti
L idea centrale dei metodi diretti Ł l idea dell eliminazione,
attribuita comunemente a Gauss, anche se pare fosse nota
precedentemente. L idea, come Ł noto, consiste nel ricavare
(eliminare) da una fissata equazione una particolare incognita e
nella sua sostituzione nelle equazioni rimanenti.
La sostituzione diminuisce la dimensione del problema da n a n-
1. Iterando il procedimento, si riduce il problema originario ad un
problema ad una dimensione, in una sola incognita. Determinata
tale incognita, le altre componenti della soluzione sono
successivamente ottenute mediante una procedura di sostituzione
all indietro.
Dal punto di vista numerico, si tratta, da una parte di organizzare
la procedura in maniera conveniente per essere implementata su
un calcolatore e dall altra di analizzare i limiti di applicazione del
metodo, cioŁ individuare per quali matrici il metodo pu essere
utilizzato e in quale forma il metodo si presenta piø stabile.
Se C e D sono due matrici qualunque, diagonali e non singolari,
un sistema dato nella forma A x = b Ł equivalente al seguente
sistema:
(C A D)( D-1)=(C b)
La trasformazione precedente corrisponde ad introdurre uno
scaling sulle variabili incognite e sui termini noti. E importante
realizzare il fatto che, benchØ matematicamente equivalenti,
differenti formulazioni possono dare, quando risolte con un
CAPITOLO 1. I METODI ITERATIVI 5
metodo numerico particolare, risultati numerici differenti. Si
pone, quindi, un problema di scelta dello scaling ottimale. Il
problema non Ł, in generale, di facile soluzione. Utili
suggerimenti per uno scaling conveniente possono, tuttavia,
venire da un esame del modello matematico da cui proviene il
sistema lineare; in generale, infatti, un buon scaling Ł il risultato
di una buona scelta di unit di misura.
CAPITOLO 1. I METODI ITERATIVI 6
1.3 Metodi iterativi
In alcune situazioni, per esempio nelle soluzioni di equazioni alle
derivate parziali di tipo ellittico con metodi alle differenze finite
o agli elementi finiti, i sistemi da risolvere sono sparsi e di
dimensioni tali (n=103:106) da rendere inutilizzabili, o quanto
meno inefficienti, i metodi diretti come il metodo delle
eliminazioni successive di Gauss, anche con i moderni
calcolatori. Infatti, mentre in questi casi la matrice iniziale ha un
numero di elementi non nulli p<<n2, il processo delle
eliminazioni successive del metodo di Gauss cambia le equazioni
del sistema ad ogni passo, in modo che la matrice dei
coefficienti pu diventare sempre meno sparsa e richiedere
quindi la memorizzazione di un numero eccessivo di elementi. I
metodi iterativi invece non alterano mai la matrice iniziale.
Partendo da un approssimazione iniziale x(0) essi definiscono una
successione di approssimazioni x(1), x(2), convergente, sotto
opportune ipotesi, alla soluzione x del sistema non singolare
A x=b.
Una nota classe di metodi iterativi pu venire definita
decomponendo dapprima la matrice A nella forma:
A = D + C
e riscrivendo poi il sistema come segue:
D x = -C x + b
Partendo da un generico vettore x(0) possiamo costruire la
CAPITOLO 1. I METODI ITERATIVI 7
successione x(1), x(2), innescando il procedimento iterativo D
x
(k+1)
= d(k) k = 0,1,2
dove
d(k) =-C x(k) + b
PoichØ il calcolo di x(k+1) comporta la soluzione del sistema
precedente, Ł chiaro che la scelta della matrice D, chiamata
matrice split , deve essere fatta in modo che det (D) ≠ 0 e,
soprattutto, il sistema risulti facilmente risolvibile.
Per esaminare la convergenza della successione x(1),
x(2), bisogna partire dalla relazione:
D x(k+1) =- C x(k) + b
da cui si ottiene
x
(k+1)
=-D-1 C x(k) + D-1 b
Ponendo poi
B = -D-1 C = I -D-1 A
la relazione precedente assume la nuova forma
x
(k+1)
=B x(k) + D-1 b (*)
Consideriamo x(k+1) e l errore assoluto e(k+1) ad esso associato:
e
(k+1)
= x- x(k+1) = (B x + D-1 b) (B x(k) + D-1 b) = B (x- x(k)) =
=B e(k) k = 0,1,2,
Possiamo allora scrivere :
e
(0)
= x - x(0)
e
(1)
= B e(0)
e
(2)
= B e(1)=BB e(0) = B2 e(0)
e
(3)
= B e(2)=BB2 e(0) = B3 e(0)
CAPITOLO 1. I METODI ITERATIVI 8
e
(k+1)
= B e(k)=BBk e(0) = Bk+1 e(0)
Per capire da cosa dipende la convergenza del processo iterativo
(*) bisogna ricordare il seguente teorema:
Sia B una matrice quadrata di ordine n. Allora :
lim Bm = O (matrice nulla)
m→∞
se e solo se ρ(B) < 1, dove ρ(B) denota il raggio spettrale di B.
Questo teorema ci permette di affermare che il processo iterativo
(*) Ł convergente, cioŁ lim e(k) = 0, se e solo se:
k→∞
ρ(I - D-1A) < 1.
PoichØ per ogni norma di matrice compatibile con una di vettore
abbiamo:
ρ(B) ≤ B
il metodo iterativo in questione Ł senz altro convergente quando
I - D-1A <1.
L idea di sdoppiare la matrice nella forma A = D + C conduce
alla costruzione di un metodo iterativo, definito dalla (*), che
comporta ad ogni passo la soluzione del sistema non singolare
D x(k+1) = d(k) k = 0,1,2
Ovviamente, affinchØ questo procedimento possa presentare dei
vantaggi il sistema:
D x(k+1) = d(k) k = 0,1,2
deve risultare piø semplice del sistema di partenza A x = b
CAPITOLO 1. I METODI ITERATIVI 9
1.4 Metodo di Jacobi
Sia D la matrice split :
D=
supponendo a11
0, i = 0,1,2, .,n. Se qualche a11 risulta nullo,
allora prima di procedere allo sdoppiamento dobbiamo permutare
le righe di A i modo che i nuovi elementi della diagonale siano
tutti nulli. Quando A Ł non singolare ci Ł sempre possibile.
Ovviamente tali scambi devono essere fatti anche sui
corrispondenti elementi di b.
Scriviamo il sistema A x = b esplicitamente:
a11 x1+ a12 x2+ + a1n xn=b1
a21 x1+ a22 x2+ + a2n xn=b2
an1 x1+ an2 x2+ + ann xn=bn
con le equazioni ordinate in modo che a11 0, i = 1,2, .,n. Tale
sistema pu anche essere riscritto nella forma:
a 11 0 0 0 .0
0 a22 0 0..... 0
0 0 a33 0. ....0
..
0 0 0 0 ...ann
CAPITOLO 1. I METODI ITERATIVI 10
=
i
x
ii
n
jj
jiji
a
xab
≠=
−
11
i = 1,2, .,n.
Il metodo di Jacobi consiste nel calcolare, nota
un approssimazione iniziale x(0) ( oppure prendendo x(0) = 0), le
approssimazioni successive ( )Tknkkk xxxx )1()1(2)1(1)1( ,....,, ++++ = k =
0,1,2 . per mezzo della relazione:
Una condizione sufficiente per la convergenza Ł:
≠=
>
n
ijj
ijij
aa
,1
i=1, .,n
Ricordando la relazione (*) del paragrafo 1.3 possiamo affermare
che il metodo risulter convergente se e solo se
ρ(BJ) < 1 dove BJ = (I - D-1 A)
e ρ(BJ) Ł il raggio spettrale della matrice di iterazione (I - D-1 A).
=+ )( k
i
x
ii
n
jj
k
jiji
a
xab
≠=
−
11
)(