IntroduzioneUn sistema di equazioni algebriche lineari ha la forma
seguente:
le N incognite x1,x2,...,xN sono legate da M equazioni linearile quantità aij e bi (i=1,2,...,M e j=1,2,...,N) sono notese N=M il numero di equazioni è pari al numero delle
incognite ed esiste la possibilità che il sistema abbia un’unica soluzione
nel caso M=N la soluzione non è unica se:una equazione è una combinazione lineare delle altre
(degenerazione di riga) tutte le equazioni contengono certe variabili nelle stesse
combinazioni lineari (degenerazione di colonna)
MNMNMMM
NN
NN
bxaxaxaxa
bxaxaxaxa
bxaxaxaxa
332211
22323222121
11313212111
Alcuni esempiDegenerazione di riga:
il primo sistema è impossibile, il secondo è indeterminato
Degenerazione di colonna:
le variabili x2 e x3 compaiono sempre nella combinazione lineare x2-x3
il sistema in esame è impossibile
0422
15
32
321
321
321
xxx
xxx
xxx
933
125
3
21
321
21
xx
xxx
xx
3227
1
52
321
321
321
xxx
xxx
xxx
Soluzioni numericheNella soluzione numerica di un sistema si
possono presentare due tipi di problemi:gli errori di arrotondamento del computer
possono rendere alcune equazioni (che inizialmente non lo sono!) linearmente dipendenti, rendendo il sistema irrisolvibile
gli errori di arrotondamento che si accumulano nel corso della risoluzione numerica possono avere effetti significativi sulla soluzione, che risulterà errata è sempre opportuno verificare la correttezza della
soluzione andando a sostituirla nelle equazioni di partenza
questa situazione tende a verificarsi se N è grande e se il sistema è vicino ad essere singolare
MatriciUn sistema lineare può scriversi nella forma:
dove:
La matrice A è detta “matrice dei coefficienti”il primo indice di ciascun elemento aij è quello
di riga, il secondo indice è quello di colonnaIl vettore b è detto “vettore dei termini noti”
bxA
MNMM
N
N
aaa
aaa
aaa
A
21
22221
11211
Nx
x
x
x
2
1
Mb
b
b
b
2
1
Definizioni (1)Matrice quadrata: M=NSia A una matrice quadrata:
diagonale principale: è formata dagli elementi con indici uguali a11,a22,...,aNN
diagonale secondaria: è formata dagli elementi con indici la cui somma è N+1: a1N,a2N-1,...,aN1
Matrice diagonale: tutti gli elementi al di fuori della diagonale principale sono nulli
Matrice simmetrica: aik=aki
NNNN
N
N
aaa
aaa
aaa
A
21
22221
11211
diagonale principale
diagonale secondaria
Definizioni (2)Matrice unità:
è una matrice quadrata che ha uguali a 1 tutti gli elementi della diagonale principale e uguali a 0 tutti gli altri elementi
Simbolo di Kronecker:
gli elementi della matrice unità possono essere rappresentati con il simbolo di Kronecker
100
010
001
I
ji
jiij
se 0
se 1
Definizioni (3)Matrice trasposta:
data una matrice A, la trasposta AT è la matrice che si ottiene da A scambiando le righe con le colonne
se A è una matrice MN, AT sarà una matrice NM
se A è quadrata, anche AT è quadrata la matrice trasposta della trasposta è nuovamente la
matrice di partenza: (AT)T=Auna matrice quadrata simmetrica coincide con la sua
trasposta
MNNN
M
M
T
MNMM
N
N
aaa
aaa
aaa
A
aaa
aaa
aaa
A
21
22212
12111
21
22221
11211
Somma e differenza tra matriciLa somma tra due matrici A e B si può effettuare solo se le
due matrici hanno le stesse dimensioni la matrice somma ha le stesse dimensioni di A e di B
Detta C=A+B la matrice somma, i suoi elementi si ottengono sommando quelli di A e B:cik=aik+bik
Data una matrice A, si definisce la matrice opposta –A tale che A+(-A)=0 la matrice nulla (0) è una matrice i cui elementi sono tutti
nulli la matrice opposta ha le stesse dimensioni della matrice A gli elementi della matrice opposta sono gli opposti di quelli
della matrice A: (-a)ik=-aik Analogamente alla somma, la differenza tra due matrici A
e B si può calcolare solo se A e B hanno le stesse dimensioni la matrice differenza ha le stesse dimensioni di A e B e si
ottiene sommando la matrice A con l’opposta di BDetta D=A-B la matrice differenza, i suoi elementi sono:
dik=aik-bik
Prodotto di una matrice per un numeroData una matrice A ed un numero r, il
prodotto del numero r per la matrice A è una matrice delle stesse dimensioni di A, i cui elementi sono pari agli elementi di A moltiplicati per il numero r:
MNMM
N
N
MNMM
N
N
rarara
rarara
rarara
ArrA
aaa
aaa
aaa
A
21
22221
11211
21
22221
11211
Prodotto tra matrici (1)Il prodotto tra matrici si può effettuare solo
se il numero di colonne della prima matrice è uguale al numero di righe della seconda matrice
Data una matrice A di tipo MP ed una matrice B di tipo PN, la matrice prodotto C=AB sarà una matrice di tipo MNsi noti che in generale AB ≠ BA (sono due
matrici di tipo diverso, ammesso che esistano entrambe!)
Gli elementi di C=AB si calcolano con la regola del prodotto righe per colonne:
P
jjkijik bac
1
Prodotto tra matrici (2)Regola del prodotto:
MNMM
ik
N
PNPkPP
Nk
Nk
MPMM
iPii
P
P
ccc
c
cc
ccc
bbbb
bbbb
bbbb
aaa
aaa
aaa
aaa
21
2221
11211
21
222221
111211
21
21
22221
11211
PkiPkiki
P
jjkijik babababac
22111
riga i
colonna k
Determinante (1)Il determinante può essere definito solamente
nel caso di matrici quadratePer una matrice quadrata 11 (del primo
ordine) il determinante è pari all’unico elemento della matrice:
Per matrici quadrate di ordine superiore il determinante si definisce in maniera induttiva, partendo dalla definizione data nel caso della matrice di ordine 1
1111 det aAAaA
Determinante (2)A seconda della sua posizione nella matrice A, un elemento aik di
una matrice può essere:di classe pari (dispari) se i+k è pari (dispari)
A ciascun elemento aik viene associata una matrice quadrata, che si ottiene dalla matrice A eliminando la i-esima riga e la k-esima colonna il determinante di tale matrice, Mik, si chiama minore
complementare dell’elemento aik
si chiama complemento algebrico Cik dell’elemento aik il minore complementare Mik moltiplicato per +1 o -1 a seconda che aik sia di classe pari o dispari:
NNNkNN
iNikii
iNk
Nk
aaaa
aaaa
aaaa
aaaa
A
21
21
22221
111211
ikki
ik MC 1
Determinante (3)Il determinante di una matrice quadrata è
definito come somma dei prodotti degli elementi di una linea (riga o colonna) per i corrispondenti complementi algebrici
Il valore del determinante non dipende dalla linea scelta
Una matrice con determinante nullo è detta singolare
)arbitrario ( det
)arbitrario ( det
1
1
iCaA
kCaA
N
kikik
N
iikik
Determinante di una matrice 22Calcoliamo il determinante di una matrice 22:
Sviluppando il determinante lungo la prima riga si ha:
Il determinante della matrice 22 si calcola come il prodotto degli elementi della diagonale principale meno il prodotto degli elementi della diagonale secondaria
2221
1211
aa
aaA
21122211213
12222
1112121111 11det aaaaaaaaCaCaA
Determinante di una matrice 33 Calcoliamo il determinante di una matrice 33:
Sviluppando il determinante lungo la prima riga si ha:
Regola empirica: il determinante di una matrice 33 si calcola sommando i prodotti degli elementi della diagonale principale e delle sue parallele e sottraendo i prodotti degli elementi della diagonale secondaria e delle sue parallele
333231
232221
131211
aaa
aaa
aaa
A
312213332112322311322113312312332211
312232211331233321123223332211
3231
2221413
3331
2321312
3332
2322211
131312121111
111
det
aaaaaaaaaaaaaaaaaa
aaaaaaaaaaaaaaa
aa
aaa
aa
aaa
aa
aaa
CaCaCaA
Alcune proprietà dei determinantiIl determinante della matrice trasposta è uguale al
determinante della matrice di partenza: detAT=detASe tutti gli elementi di una riga (o di una colonna) sono
nulli, il determinante è nulloScambiando tra loro due righe (o due colonne) il
determinante cambia segnoSe due righe (o due colonne) parallele sono
proporzionali, il determinante è nulloSe gli elementi di una riga (o di una colonna) vengono
moltiplicati per un numero r, il determinante viene anch’esso moltiplicato per r
Il determinante di una matrice quadrata non cambia quando agli elementi di una linea si aggiungono i corrispondenti elementi di una linea parallela moltiplicati per un numero arbitrario r
Date due matrici quadrate A e B, se C=AB, allora risulta anche detC=detAdetB
Matrici triangolariUna matrice quadrata si dice triangolare se tutti i
suoi elementi al di sopra o al di sotto della diagonale principale sono nullimatrice triangolare superiore (alta) se aik=0 con i>k
matrice triangolare inferiore (bassa) se aik=0 con i<k
Il determinante di una matrice triangolare è pari al prodotto degli elementi della diagonale principale:
NNNNNN
N
N
U
aaa
aa
a
T
a
aa
aaa
T
21
2221
11
L
222
11211
0
00
00
0
N
iiiNN aaaaT
12211det
Matrice inversaData una matrice quadrata A di ordine N, si dice
inversa (se esiste) la matrice A-1 tale che:
dove I è la matrice unità di ordine NOgni matrice che ammette un’inversa si dice
invertibileSe la matrice inversa esiste è anch’essa di ordine N
l’inversa della matrice unità è ancora la matrice unitàTeorema: se la matrice inversa esiste, essa è unicaTeorema: ogni matrice non singolare possiede una
matrice inversase una matrice ha determinante nullo, allora non è
invertibile
AAIAA 11
Metodo di Gauss-Jordan (1)Il problema della risoluzione di un sistema
lineare è strettamente connesso a quello dell’inversione della matrice dei coefficientii due problemi possono essere affrontati insieme
Sia quindi data una matrice A di NN coefficienti ed un vettore b di N termini noti
Vogliamo trovare un vettore x ed una matrice Y tali che:
il vettore x è la soluzione di un sistema lineare di N equazioni in N incognite
la matrice Y è l’inversa della matrice A
IAY
bAx
Metodo di Gauss-Jordan (2)Osservazioni preliminari:
se si scambiano tra loro due righe di A, e le corrispondenti due righe di b e di I (che in seguito allo scambio non sarà più la matrice unità), le soluzioni x e Y non cambiano
se sostituiamo una riga di A con una combinazione lineare della riga stessa e di una qualsiasi altra riga, e facciamo la stessa cosa nelle righe di b e di I (che, ancora una volta, in seguito allo scambio non sarà più la matrice unità), le soluzioni x e Y non cambiano
se si scambiano due colonne di A, si ottengono le stesse soluzioni solo se simultaneamente vengono scambiate le righe corrispondenti di x e di Y lo scambio di colonne comporta uno scambio delle righe nella
soluzioneIl metodo di eliminazione di Gauss-Jordan sfrutta
tutte queste proprietà
Un esempio (1)Consideriamo il sistema:
La soluzione è il vettore x1=3, x2=2, x3=1Scambiamo la prima e la seconda riga di A e la prima e la
seconda riga di b
è evidente che il sistema è rimasto lo stesso, e quindi anche la soluzione non è cambiata
83
52
342
321
321
321
xxx
xxx
xxx
8
5
3
b
311
112
421
A
83
342
52
321
321
321
xxx
xxx
xxx
8
3
5
b
311
421
112
A
Un esempio (2)Sostituiamo ora la prima riga di A e di b con
la somma della prima riga e della seconda riga
Ancora una volta la soluzione del sistema non cambia (x1=3, x2=2, x3=1)
8
5
3
b
311
112
421
A
8
5
8
b
311
112
313
A
83
52
833
321
321
321
xxx
xxx
xxx
Un esempio (3)Infine scambiamo la seconda e la terza colonna di A:
Adesso la soluzione è diventata x1=3, x2=1, x3=2, cioè si sono scambiati x2 con x3
Per ripristinare la soluzione precedente è necessario scambiare la seconda e la terza riga del vettore x
8
5
3
b
311
112
421
A
8
5
3
b
131
112
241
A
83
52
324
321
321
321
xxx
xxx
xxx
Metodo di Gauss-Jordan (3)In ciascuna iterazione la matrice A, il vettore b e la
matrice I vengono modificati sfruttando le proprietà enunciate in precedenza
Alla n-esima iterazione le equazioni da risolvere avranno la forma:
se si riesce a trasformare la matrice A nella matrice unità, ossia se An=I, allora sarà:
cioè il vettore dei termini noti modificato e la matrice unità modificata forniranno la soluzione del problema
nn
nn
IYA
bxA
n
n
IY
bx
Esempio (1)Vediamo come si applica il metodo di eliminazione
di Gauss-Jordan nell’esempio studiato in precedenza:
Prima iterazione: dividiamo la prima riga per l’elemento a11 (pivot) in modo che il primo elemento della diagonale principale sia 1in questo caso particolare a11=1 e le matrici restano
invariate
100
010
001
8
5
3
b
311
112
421
IA
100
010
001
8
5
3
b
311
112
421
111 IA
Esempio (2)Seconda iterazione:
lo scopo è quello di far sì che la prima colonna riproduca la prima colonna della matrice unità seconda riga seconda riga – a21prima riga terza riga terza riga – a31prima riga
100
010
001
8
5
3
b
311
112
421
111 IA
101
012
001
5
1
3
b
710
950
421
222 IA
Esempio (3)Terza iterazione:
lo scopo è quello di far sì che anche il secondo elemento della diagonale principale di A sia pari a 1
dividiamo quindi la seconda riga per l’elemento pivot a22
101
012
001
5
1
3
b
710
950
421
222 IA
101
05/15/2
001
5
5/1
3
b
710
5/910
421
333 IA
101
05/15/2
001
5
5/1
3
b
710
5/910
421
333 IA
Esempio (4)Quarta iterazione:
lo scopo è quello di far sì che la seconda colonna riproduca la seconda colonna della matrice unità prima riga prima riga – a12seconda riga terza riga terza riga – a32seconda riga
15/15/3
05/15/2
05/25/1
5/26
5/1
5/13
b
5/2600
5/910
5/201
444 IA
15/15/3
05/15/2
05/25/1
5/26
5/1
5/13
b
5/2600
5/910
5/201
444 IA
Esempio (5)Quinta iterazione:
lo scopo è quello di far sì che anche il terzo elemento della diagonale principale di A sia pari a 1
dividiamo quindi la terza riga per l’elemento pivot a33
26/526/126/3
05/15/2
05/25/1
1
5/1
5/13
b
100
5/910
5/201
555 IA
26/526/126/3
05/15/2
05/25/1
1
5/1
5/13
b
100
5/910
5/201
555 IA
Esempio (6)Sesta iterazione:
lo scopo è quello di far sì che la terza colonna riproduca la seconda colonna della matrice unità prima riga prima riga – a13terza riga seconda riga seconda riga – a23terza riga
26/526/126/3
26/926/726/5
13/113/513/2
1
2
3
b
100
010
001
666 IA
Alcune considerazioniNell’esempio precedente il vettore b6 rappresenta la
soluzione x del sistema Ax=b e la matrice I6 rappresenta la matrice Y=A-1
Nell’applicare il metodo di Gauss-Jordan potrebbero nascere dei problemi se un elemento pivot dovesse essere nullo in tal caso gli elementi della riga corrispondente verrebbero
divisi per zeroPer evitare situazioni in cui un elemento pivot è nullo, si
sfruttano le proprietà legate alla permutazione di righe e colonne delle matrici, in modo che nella posizione di pivot si venga a trovare un elemento di matrice “opportuno”“partial pivoting”: si sfruttano solo le proprietà legate alla
permutazione delle righe“full pivoting”: si sfruttano le proprietà legate sia alla
permutazione delle righe che delle colonneuna volta trovate le soluzioni del problema, occorre poi
effettuare le opportune permutazioni inverseuna buona strategia è quella di scegliere volta per volta come
elemento pivot quello più grande in valore assoluto, effettuando le permutazioni necessarie per portarlo nella posizione di pivot
Metodo di eliminazione di Gauss Il metodo di eliminazione di Gauss trasforma una matrice A
in una matrice triangolare Una volta effettuata questa operazione, il sistema Ax=b
può essere risolto per sostituzioneAlla fine della procedura, si ha:
Il sistema così ottenuto può essere risolto per sostituzione:
NNNN
NN
NN
NNN
N
N
bxa
bxaxa
bxaxaxa
bAx
b
b
b
b
a
aa
aaa
A
22222
11212111
2
1
222
11211
00
0
N
ijjiji
iii
NN
NN
xaba
x
a
bx
1
1
Esempio (1)Consideriamo ancora il sistema precedente:
Prima iterazione:facciamo in modo che a21=0 e a31=0
seconda riga seconda riga – (a21/a11)prima riga
terza riga terza riga – (a31/a11)prima riga
83
52
342
8
5
3
b
311
112
421
321
321
321
xxx
xxx
xxx
A
57
195
342
5
1
3
b
710
950
421
32
32
321
11
xx
xx
xxx
A
Esempio (2)Seconda iterazione:
facciamo in modo che a32=0
terza riga terza riga – (a32/a22)seconda riga
57
195
342
5
1
3
b
710
950
421
32
32
321
11
xx
xx
xxx
A
5
26
5
26195
342
5/26
1
3
b
5/2600
950
421
3
32
321
22
x
xx
xxx
A
Esempio (3)Il sistema può quindi essere risolto per
sostituzione partendo dall’ultima equazione:
3423
25
19
1
321
32
3
xxx
xx
x
Considerazioni sul metodo di GaussAnche in questo caso bisogna evitare le
situazioni in cui l’elemento pivot è nullo, sfruttando le proprietà delle matrici
Il metodo di eliminazione di Gauss è più veloce rispetto a quello di Gauss-Jordannegli algoritmi in cui viene applicato il metodo di
Gauss-Jordan sono richieste N3 operazioni per ciclo
negli algoritmi in cui si applica il metodo di Gauss sono richieste N3/3 operazioni per ciclo
Il metodo di eliminazione di Gauss non fornisce la matrice inversa A-1
Tale matrice può essere calcolata risolvendo gli N sistemi Axi=bi dove bi è il vettore contenente l’i-esima colonna della matrice unità
Decomposizione LU (1)Supponiamo che sia possibile scrivere la
matrice A come prodotto A=LU dove:L è una matrice triangolare inferioreU è una matrice triangolare superiore
Se tale decomposizione è possibile, il sistema Ax=b si trasforma nel modo seguente:
NN
N
N
NNNNNNNNNN
N
N
aaaa
aaaa
aaaa
LUA
000
000
000
22322
1131211
321
2221
11
321
2232221
1131211
bUxLbxLUbAx
Decomposizione LU (2)Per risolvere il sistema di partenza occorre dunque
risolvere i due sistemi:Ly = bUx = y
Il vantaggio di questo approccio è che entrambi i sistemi hanno una matrice dei coefficienti di forma triangolare, e possono essere risolti per sostituzione:
Niyby
by
i
jjiji
iii ,,3,2
1 1
1
11
11
1,,2,1 1
1
NNixyx
yx
N
ijjiji
iii
NN
NN
Decomposizione LU (3)Se vale la decomposizione A=LU si ha, dalla
regola del prodotto:
in totale si possono scrivere N2 equazioni di questo tipociascuna equazione contiene un numero di termini
che dipende da i e j (km=0 se k<m e km=0 se k>m)
le incognite sono N2+N (le ij e le ij non nulle)in ciascuna matrice L e U ci sono N(N+1)/2
incogniteper effettuare la decomposizione occorre fissare N
incognite
NjiNji
N
kjikjikij aaaa
221
11
Algoritmo di Crout (1)Si fissano i termini diagonali della matrice L
ponendo:
Partiamo dalle equazioni che legano i coefficienti ij, ij e aij ed eliminiamo i termini nulli tenendo conto che km=0 se k<m e km=0 se k>m:
Nell’algoritmo di Crout si sfruttano queste relazioni per ricavare, in maniera iterativa, i coefficienti ij e ij
Niii ,,2,1 1
N
k
jjijjijiij
jjiijijiij
ijiijijiij
kjikij
jia
jia
jia
a1
2211
2211
2211
Algoritmo di Crout (2)Fissiamo l’indice di colonna j=1,2,...,N.
Facciamo variare l’indice di riga i nell’intervallo i=1,2,...,j:
Facciamo ora variare l’indice di riga i nell’intervallo i=j+1,j+2,...,N:
I coefficienti che compaiono nelle equazioni sono stati determinati nelle iterazioni precedenti
nulla siconsiderar da è sommatoria la 1 se 1
1
1
12211
ia
a
i
kkjikijij
ij
i
kkjikijiijijiij
nulla) siconsiderar da è sommatoria la 1 (se 1 1
1
1
12211
ja
a
j
kkjikij
jjij
jjij
j
kkjikjjijjijiij
Esempio (1)Applichiamo l’algoritmo di Crout per
decomporre la matrice:
Alla fine del procedimento determineremo due matrici L e U tali che:
311
112
421
A
33
2322
131211
3231
21
00
0
1
01
001
ULLUA
Esempio (2)L’algoritmo di Crout viene applicato colonna
per colonnasi determinano prima i valori di ij e poi quelli
di ij
Colonna 1 (j=1):
Colonna 2 (j=2):
111
11
221
11
1
3111
31
2111
21
1111
a
a
a
5
1211
5
11
5221
2
12313222
32
12212222
1212
a
a
a
Esempio (3)Colonna 3 (j=3):
Riassumendo, si ha:
e si può verificare che A=LU
5
269
5
1433
9421
4
233213313333
13212323
1313
a
a
a
5/2600
950
421
15/11
012
001
UL
Commenti sull’algoritmo di CroutAnche nell’algoritmo di Crout compaiono delle
divisioni per gli elementi pivot jj
la scelta dell’elemento pivot può essere effettuata mediante scambi di righe (“partial pivoting”)una volta calcolati tutti i ij della j-esima colonna, si
effettua la permutazione delle righe che porta sulla diagonale l’elemento voluto (di solito si sceglie il più grande)
in questo caso la matrice che viene decomposta non è la matrice originaria A, ma una matrice ottenuta permutando le righe di A: occorre quindi ricordare le permutazioni effettuate
La soluzione di un sistema lineare tramite la decomposizione LU richiede circa 1/3 delle operazioni richieste dal metodo di Gauss-Jordan
Applicazioni della decomposizione LUCalcolo della matrice inversa:
come nel caso del metodo di Gauss, la matrice inversa può essere calcolata risolvendo gli N sistemi Axi=bi dove bi è il vettore contenente l’i-esima colonna della matrice unità
Calcolo del determinante:si sfruttano le proprietà del determinante di un
prodotto tra matrici e del determinante delle matrici triangolari:
N
iii
N
iii
N
iiiULA
111
detdetdet
Sistemi tridiagonali (1)Si tratta di sistemi in cui la matrice dei
coefficienti ha elementi non nulli solo sulla diagonale principale e nelle posizioni adiacenti ad essa:
In generale gli elementi sulla diagonale principale sono più grandi rispetto a quelli fuori diagonale:
NN
A
000
000
00
000
33
222
11
iii
Sistemi tridiagonali (2)La soluzione di un sistema tridiagonale richiede un
numero di operazioni minore rispetto a quella di un sistema normale tipicamente per risolvere un sistema normale sono
richieste O(N3) operazioni, mentre per risolvere un sistema tridiagonale ne occorrono O(N)
Un algoritmo per la soluzione di un sistema tridiagonale richiede meno memoria rispetto ad uno stesso algoritmo applicato ad un sistema normalenon è necessario immagazzinare le informazioni complete
sulla matrice dei coefficienti, ma solo quelle relative ai coefficienti non nulli
Si può sviluppare un algoritmo per la soluzione di sistemi tridiagonali che sfrutta il metodo eliminazione di Gauss la soluzione si ricava per sostituzionepoiché solitamente gli elementi sulla diagonale sono diversi
da zero e più grandi in modulo di quelli fuori diagonale, in genere non ci sono problemi di scelta dell’elemento pivot
Esempio (1)Consideriamo il seguente sistema
tridiagonale:
la cui soluzione è x1=2, x2=3, x3=1, x4=1Risolviamo il sistema trasformando la matrice
A in una matrice triangolare con il metodo di eliminazione di Gauss
13
8
15
7
10300
3810
0251
0015
13103
838
1525
75
43
432
321
21
BA
xx
xxx
xxx
xx
Esempio (2)Prima iterazione:
seconda riga seconda riga - 2/1prima riga
13
8
15
7
10300
3810
0251
0015
BA
13
8
5/68
7
10300
3810
025/260
0015
11 BA
13
8
5/68
7
10300
3810
025/260
0015
11 BA
Esempio (3)Seconda iterazione:
terza riga terza riga - 3/2seconda riga
13
13/70
5/68
7
10300
313/10900
025/260
0015
22 BA
13
13/70
5/68
7
10300
313/10900
025/260
0015
22 BA
Esempio (4)Terza iterazione:
quarta riga quarta riga - 4/3terza riga
109/1207
13/70
5/68
7
109/1207000
313/10900
025/260
0015
33 BA
Esempio (5)A questo punto il sistema può essere risolto
per sostituzione:
In pratica, data la struttura tridiagonale, in ogni iterazione viene cambiato soltanto un elemento della matrice A
275
175
325
68
26
5
5
682
5
26
1313
70
109
13
13
703
13
109
1109
1207
109
1207
2121
3232
4343
44
xxxx
xxxx
xxxx
xx
Metodi iterativiUn sistema di equazioni lineari può essere
risolto in maniera iterativaSi parte da una soluzione «di prova» e,
attraverso iterazioni successive, si arriva alla soluzione del sistema
Non sempre le procedure iterative convergono verso la soluzione del sistemaSe la matrice dei coefficienti soddisfa una serie
di condizioni opportune, la convergenza è garantita
Una volta stabilita la convergenza della procedura iterativa, occorre stabilire un criterio per arrestare le iterazioni
Metodo di Jacobi (1)Consideriamo come esempio il sistema
seguente:
la cui soluzione è x=2, y=4, z=3Le equazioni del sistema possono essere
riscritte nel modo seguente:
1552
2184
74
zyx
zyx
zyx
5
2158
4214
7
yxz
zxy
zyx
Metodo di Jacobi (2)Si può quindi pensare, data la soluzione (xk,yk,zk)
dopo k iterazioni, di valutare la soluzione (xk+1,yk+1,zk+1) all’iterazione successiva come:
Partendo da una soluzione di prova (x0,y0,z0) ci si aspetta che la procedura iterativa converga verso la soluzione esatta del sistema in un numero ragionevole di passi
5
2158
4214
7
1
1
1
kkk
kkk
kkk
yxz
zxy
zyx
Metodo di Jacobi (3)Partendo dai valori (x0,y0,z0)=(1,2,2) il metodo di Jacobi
converge verso la soluzione in 10 iterazioni (con una precisione alla quinta cifra significativa)
Iterazione x y z
0 1 2 2
1 1.75 3.375 3
2 1.8438 3.875 3.025
3 1.9625 3.925 2.9625
4 1.9906 3.9766 3
5 1.9941 3.9953 3.0009
6 1.9986 3.9972 2.9986
7 1.9996 3.9991 3
8 1.9998 3.9998 3
9 1.9999 3.9999 2.9999
10 2 4 3
Metodo di Jacobi (4)Riscriviamo il sistema di partenza
riordinando le equazioni nel modo seguente:
Ricavando x dalla prima equazione, y dalla seconda e z dall’ultima possiamo scrivere le formule di Jacobi nel modo seguente:
74
2184
1552
zyx
zyx
zyx
yxz
zxy
zyx
478
4212
515
kkk
kkk
kkk
yxz
zxy
zyx
478
4212
515
1
1
1
Metodo di Jacobi (5)Partendo dai valori (x0,y0,z0)=(1,2,2) e utilizzando le
formule ricorsive ricavate in precedenza, il metodo di Jacobi stavolta non converge verso la soluzione!
Iterazione x y z
0 1 2 2
1 -1.5 3.375 5
2 6.6875 2.5 16.375
3 34.6875 8.0156 -17.25
4 -46.6172 17.8125 -123.7344
5 -307.9227 -36.1504 211.2813
6 502.6279 -124.9297 1202.5684
Metodo di Gauss-Seidel (1)Consideriamo lo stesso sistema dell’esempio
precedente
L’idea alla base del metodo di Gauss-Seidel è quella di velocizzare la ricerca della soluzione sfruttando i valori di ciascuna incognita, man mano che questi vengono calcolati
1552
2184
74
zyx
zyx
zyx
5
2158
4214
7
yxz
zxy
zyx
Metodo di Gauss-Seidel (2)La formula che permette di passare dalla k-
esima alla (k+1)-esima iterazione sarà dunque la seguente:
Anche in questo caso si parte da una soluzione di prova (x0,y0,z0) e si procede iterativamente
5
2158
4214
7
111
11
1
kkk
kkk
kkk
yxz
zxy
zyx
Metodo di Gauss-Seidel (3)Partendo dai valori (x0,y0,z0)=(1,2,2) e utilizzando
le formule ricorsive ricavate in precedenza, il metodo di Gauss-Seidel converge verso la soluzione in 6 iterazioni (con una precisione alla quinta cifra significativa)
Iterazione x y z
0 1 2 2
1 1.75 3.75 2.95
2 1.95 3.9688 2.9862
3 1.9956 3.9961 2.999
4 1.9993 3.9995 2.9998
5 1.9999 3.9999 3
6 2 4 3
Matrici a diagonale dominanteUna matrice A è detta a diagonale dominante (per
righe) in senso stretto se è verificata la condizione:
Questo significa che in ciascuna riga l’elemento che si trova sulla diagonale è in valore assoluto maggiore della somma dei valori assoluti degli altri elementiIn maniera analoga si definiscono le matrici a
diagonale dominante per colonneIn questo caso, considerata ciascuna colonna,
l’elemento che si trova sulla diagonale deve essere in valore assoluto maggiore della somma dei valori assoluti degli altri elementi
Teorema: una matrice a diagonale dominante è non singolare
NkaaN
kjj
kjkk 1 1
EsempioConsideriamo il sistema visto negli esempi
precedenti:
Come si può facilmente verificare, la matrice A è strettamente a diagonale dominante (per righe)
Se però riscriviamo il sistema scambiando la prima e l’ultima equazione si ha:
In questo caso, invece, la matrice A non è più strettamente a diagonale dominante
1552
2184
74
zyx
zyx
zyx
512
184
114
A
74
2184
1552
zyx
zyx
zyx
114
184
512
A
Condizione sufficiente di convergenzaSe la matrice dei coefficienti A è strettamente a
diagonale dominante, allora il sistema Ax=b ammette un’unica soluzione
Si può dimostrare che in tal caso la procedura di Jacobi converge verso la soluzione, comunque si scelga il vettore x0 di partenzaLo stesso risultato vale anche per il metodo di
Gauss-SeidelIn generale, il metodo di Gauss-Seidel converge
più velocemente rispetto al metodo di Jacobi, e quindi viene preferito
Vi sono però dei casi in cui il metodo di Jacobi converge, mentre il metodo di Gauss-Seidel non converge!
Formalismo generaleDato il sistema Ax=b, le formule generali
utilizzate per le iterazioni di Jacobi e di Gauss-Seidel sono le seguenti:
Iterazione di Jacobi:
Iterazione di Gauss-Seidel:
jj
kNjN
kjjj
kjjj
kjjk
j a
xaxaxaxabx
)()(11
)(11
)(111
jj
N
jii
kijij
kj a
xab
x
1
)(
1
jj
kNjN
kjjj
kjjj
kjjk
j a
xaxaxaxabx
)()(11
)1(11
)1(111
jj
j
i
N
ji
kiji
kijij
kj a
xaxab
x
1
1 1
)()1(
1
Criterio di convergenzaPer fermare la procedura iterativa occorre che
il vettore x(k+1) sia sufficientemente vicino al vettore x(k)
Si va quindi a valutare la distanza tra i due vettori e si richiede che essa sia minore di un valore di soglia prefissato ε
Possibili condizioni di convergenza:Condizione sulla norma euclidea:
Condizione sulla norma generalizzata:
La condizione sulla norma generalizzata è più forte di quella sulla norma euclidea
N
j
kj
kj
kk xxxx1
2)()1()1(
N
j
kj
kj
kk xxxx1
)()1(
1
)1(
Top Related