Introduzione
 Un sistema di equazioni algebriche lineari ha la forma seguente:
a11 x1  a12 x2  a13 x3    a1N xN  b1
a21x1  a22 x2  a23 x3    a2 N x N  b2

aM 1 x1  aM 2 x2  aM 3 x3    aMN xN  bM
 le N incognite x1,x2,...,xN sono legate da M equazioni lineari
 le quantità aij e bi (i=1,2,...,M e j=1,2,...,N) sono note
 se 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)
Alcuni esempi
 Degenerazione di riga:
 x1  x2  2 x3  3

 x1  x2  5 x3  1

2 x1  2 x2  4 x3  0
 x1  x2  3

 x1  5 x2  2 x3  1

 3x1  3x2  9
 il primo sistema è impossibile, il secondo è indeterminato
 Degenerazione di colonna:
 2 x1  x2  x3  5

 x1  x2  x3  1

7 x1  2 x2  2 x3  3
 le variabili x2 e x3 compaiono sempre nella combinazione
lineare x2-x3
 il sistema in esame è impossibile
Soluzioni numeriche
 Nella 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
Matrici
 Un sistema lineare può scriversi nella forma:
dove:
 a11 a12  a1N 


 a21 a22  a2 N 

A





aM 1 aM 2  aMN 
A x  b
 x1 
 
 x2 
x 
 
 
 x N 
 b1 
 
 b2 
b 
 
 
bM 
 La matrice A è detta “matrice dei coefficienti”
 il primo indice di ciascun elemento aij è quello di riga, il
secondo indice è quello di colonna
 Il vettore b è detto “vettore dei termini noti”
Definizioni (1)
 Matrice quadrata: M=N
 Sia A una matrice quadrata:
 a11 a12  a1N 


 a21 a22  a2 N 

A





a N 1 a N 2  a NN 
diagonale principale
diagonale secondaria
 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
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
1 0  0


0 1  0 

I 
 



0 0  1
 Simbolo di Kronecker:
1 se i  j
 ij  
0 se i  j
 gli elementi della matrice unità possono essere rappresentati
con il simbolo di Kronecker
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 MN, AT sarà una matrice NM
 a11 a12  a1N 
 a11



 a21 a22  a2 N 
 a12
T
 A 
A







aM 1 aM 2  aMN 
a1N
a21
a22

a2 N
 aM 1 

 aM 2 



 aMN 
 se A è quadrata, anche AT è quadrata
 la matrice trasposta della trasposta è nuovamente la matrice di
partenza: (AT)T=A
 una matrice quadrata simmetrica coincide con la sua trasposta
Somma e differenza tra matrici
 La 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 B
 Detta D=A-B la matrice differenza, i suoi elementi sono:
 dik=aik-bik
Prodotto di una matrice per un numero
 Data 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:
 a11 a12  a1N 
 ra11 ra12  ra1N 




 a21 a22  a2 N 
 ra21 ra22  ra2 N 
  rA  Ar  

A










aM 1 aM 2  aMN 
ra M 1 ra M 2  raMN 
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 MP ed una matrice B di
tipo PN, la matrice prodotto C=AB sarà una matrice
di tipo MN
 si noti che in generale AB ≠ BA (sono due matrici di
tipo diverso, ammesso che esistano entrambe!)
 Gli elementi di C=AB si calcolano con la regola del
prodotto righe per colonne:
P
cik   aijb jk
j 1
Prodotto tra matrici (2)
 Regola del prodotto:
P
cik   aijb jk  ai1b1k  ai 2b2k    aiPbPk
j 1
colonna k
 a11

 a21

riga i  

 ai1



a
 M1
a12

a22

ai 2

aM 2 
a1P   b11 b12  b1k  b1N   c11 c12
 
 
a2 P  b21 b22  b2 k  b2 N   c21 c22
 
 
 
 


 

aiP  
 
 
 
 
 
 
bPN  cM 1 cM 2
aMP  bP1 bP 2  bPk
cik
c1N 









cMN 
Determinante (1)
 Il determinante può essere definito solamente nel caso
di matrici quadrate
 Per una matrice quadrata 11 (del primo ordine) il
determinante è pari all’unico elemento della matrice:
A  a11   det A  A  a11
 Per matrici quadrate di ordine superiore il
determinante si definisce in maniera induttiva,
partendo dalla definizione data nel caso della matrice
di ordine 1
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
 a11 a12

 a21 a22




A
 ai1 ai 2




a
 N1 aN 2
 a1k
 a2 k

aik
 a Nk
 a1N 

 aiN 



 aiN 



 a NN 
 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:
Cik   1 M ik
i k
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
N
det A   aik Cik (k arbitrario )
i 1
N
det A   aik Cik (i arbitrario )
k 1
 Il valore del determinante non dipende dalla linea scelta
 Una matrice con determinante nullo è detta singolare
Determinante di una matrice 22
 Calcoliamo il determinante di una matrice 22:
 a11 a12 
A

a21 a22 
 Sviluppando il determinante lungo la prima riga si ha:
det A  a11C11  a12C12  a11  1 a22  a12  1 a21  a11a22  a12a21
2
3
 Il determinante della matrice 22 si calcola come il
prodotto degli elementi della diagonale principale
meno il prodotto degli elementi della diagonale
secondaria
Determinante di una matrice 33
 Calcoliamo il determinante di una matrice 33:
 a11 a12

A  a21 a22

a31 a32
a13 

a23 

a33 
 Sviluppando il determinante lungo la prima riga si ha:
det A  a11C11  a12C12  a13C13 
a11  1
2
a22
a23
a32
a33
 a12  1
3
a21 a23
a31 a33
 a13  1
4
a21 a22
a31

a32
a11 a22a33  a23a32   a12 a21a33  a23a31   a13 a21a32  a22a31  
a11a22a33  a12a23a31  a13a21a32  a11a23a32  a12a21a33  a13a22a31
 Regola empirica: il determinante di una matrice 33 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
Alcune proprietà dei determinanti
 Il determinante della matrice trasposta è uguale al determinante






della matrice di partenza: detAT=detA
Se tutti gli elementi di una riga (o di una colonna) sono nulli, il
determinante è nullo
Scambiando tra loro due righe (o due colonne) il determinante
cambia segno
Se due righe (o due colonne) parallele sono proporzionali, il
determinante è nullo
Se 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=AB, allora risulta anche
detC=detAdetB
Matrici triangolari
 Una matrice quadrata si dice triangolare se tutti i suoi
elementi al di sopra o al di sotto della diagonale
principale sono nulli
 matrice triangolare superiore (alta) se aik=0 con i>k
 matrice triangolare inferiore (bassa) se aik=0 con i<k
a11 a12  a1N 


 0 a22  a2 N 

TU  





 0
0  a NN 
0  0 
 a11


 a21 a22  0 

TL  





a N 1 a N 2  a NN 
 Il determinante di una matrice triangolare è pari al
prodotto degli elementi della diagonale principale:
N
det T  a11a22  a NN   aii
i 1
Matrice inversa
 Data una matrice quadrata A di ordine N, si dice inversa (se
esiste) la matrice A-1 tale che:
AA1  I  A1 A
dove I è la matrice unità di ordine N
 Ogni matrice che ammette un’inversa si dice invertibile
 Se 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 è unica
 Teorema: ogni matrice non singolare possiede una matrice
inversa
 se una matrice ha determinante nullo, allora non è invertibile
Metodo di Gauss-Jordan (1)
 Il problema della risoluzione di un sistema lineare è
strettamente connesso a quello dell’inversione della
matrice dei coefficienti
 i due problemi possono essere affrontati insieme
 Sia quindi data una matrice A di NN coefficienti ed un
vettore b di N termini noti
 Vogliamo trovare un vettore x ed una matrice Y tali che:
Ax  b
AY  I
 il vettore x è la soluzione di un sistema lineare di N
equazioni in N incognite
 la matrice Y è l’inversa della matrice A
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 soluzione
 Il metodo di eliminazione di Gauss-Jordan sfrutta tutte queste
proprietà
Un esempio (1)
 Consideriamo il sistema:
 x1  2 x2  4 x3  3

 2 x1  x2  x3  5

 x1  x2  3x3  8
1 2  4 


A  2  1 1 


1 1
3 
3
 
b  5
 
8
 La soluzione è il vettore x1=3, x2=2, x3=1
 Scambiamo la prima e la seconda riga di A e la prima e la seconda riga
di b
2  1 1 


A  1 2  4


1 1
3 
5
 
b  3
 
8
 2 x1  x2  x3  5

 x1  2 x2  4 x3  3

 x1  x2  3x3  8
 è evidente che il sistema è rimasto lo stesso, e quindi anche la soluzione
non è cambiata
Un esempio (2)
 Sostituiamo ora la prima riga di A e di b con la somma
della prima riga e della seconda riga
1 2  4 


A  2  1 1 


1 1
3 
3
 
b  5
 
8
3 1  3


A  2  1 1 


1 1
3 
8
 
b  5
 
8
3x1  x2  3x3  8

 2 x1  x2  x3  5

 x1  x2  3x3  8
 Ancora una volta la soluzione del sistema non cambia
(x1=3, x2=2, x3=1)
Un esempio (3)
 Infine scambiamo la seconda e la terza colonna di A:
1 2  4 


A  2  1 1 


1 1
3 
3
 
b  5
 
8
1  4 2 


A  2 1  1


1 3
1 
3
 
b  5
 
8
 x1  4 x2  2 x3  3

 2 x1  x2  x3  5

 x1  3x2  x3  8
 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
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:
An x  bn
AnY  I n
 se si riesce a trasformare la matrice A nella matrice unità, ossia se
An=I, allora sarà:
x  bn
Y  In
 cioè il vettore dei termini noti modificato e la matrice unità
modificata forniranno la soluzione del problema
Esempio (1)
 Vediamo come si applica il metodo di eliminazione di
Gauss-Jordan nell’esempio studiato in precedenza:
1 2  4 


A  2  1 1 


1 1
3 
3
 
b  5
 
8
1 0 0


I  0 1 0 


0 0 1
 Prima iterazione: dividiamo la prima riga per l’elemento a11
(pivot) in modo che il primo elemento della diagonale
principale sia 1
 in questo caso particolare a11=1 e le matrici restano invariate
1 2  4 


A1  2  1 1 


1 1
3 
3
 
b1  5
 
8
1 0 0


I1  0 1 0


0 0 1
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 – a21prima riga
terza riga  terza riga – a31prima riga
1 2  4 


A1  2  1 1 


1 1
3 
1 2  4


A2  0  5 9 


0  1 7 
3
 
b1  5
 
8
3
 
b 2   1
 
 5 
1 0 0


I1  0 1 0


0 0 1
 1 0 0


I 2    2 1 0


  1 0 1
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
1 2  4


A2  0  5 9 


0  1 7 
4 
1 2


A3  0 1  9 / 5


0  1
7 
3
 
b 2   1
 
 5 
 3 
 
b3  1 / 5
 
 5 
 1 0 0


I 2    2 1 0


  1 0 1
0
0
 1


I 3   2 / 5  1 / 5 0


  1
0
1
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 – a12seconda riga
terza riga  terza riga – a32seconda riga
4 
1 2


A3  0 1  9 / 5


0  1
7 
1 0  2 / 5


A4  0 1  9 / 5


0 0 26 / 5 
 3 
 
b3  1 / 5
 
 5 
13 / 5 


b4   1/ 5 


26 / 5
0
0
 1


I 3   2 / 5  1 / 5 0


  1
0
1
2 / 5 0
 1/ 5


I 4   2 / 5  1 / 5 0


 3 / 5  1 / 5 1
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
1 0  2 / 5


A4  0 1  9 / 5


0 0 26 / 5 
1 0  2 / 5


A5  0 1  9 / 5


0 0
1 
13 / 5 


b4   1/ 5 


26 / 5
13 / 5


b5   1 / 5 


 1 
2 / 5 0
 1/ 5


I 4   2 / 5  1 / 5 0


 3 / 5  1 / 5 1
2/5
0 
 1/ 5


I5   2 / 5
1/ 5
0 


 3 / 26  1 / 26 5 / 26
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 – a13terza riga
seconda riga  seconda riga – a23terza riga
1 0  2 / 5


A5  0 1  9 / 5


0 0
1 
1 0 0


A6  0 1 0


0 0 1
13 / 5


b5   1 / 5 


 1 
 3
 
b 6   2
 
1
2/5
0 
 1/ 5


I5   2 / 5
1/ 5
0 


 3 / 26  1 / 26 5 / 26
5 / 13 1 / 13 
 2 / 13


I 6   5 / 26  7 / 26 9 / 26


 3 / 26  1 / 26 5 / 26
Alcune considerazioni
 Nell’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
zero
 Per 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 colonne
 una volta trovate le soluzioni del problema, occorre poi effettuare le
opportune permutazioni inverse
 una 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 sostituzione
 Alla fine della procedura, si ha:
a11 a12  a1N 
 b1 


 
 0 a22  a2 N 
 b2 
 b 
A


 



 
 0
bN 
0  a NN 
a11x1  a12 x2    a1N xN  b1

a22 x2    a2 N xN  b2
Ax  b  


a NN xN  bN
 Il sistema così ottenuto può essere risolto per sostituzione:
xN 
bN
a NN
1
xi 
aii
N


b

a
x
 i  ij j 
j i 1


Esempio (1)
 Consideriamo ancora il sistema precedente:
1 2  4 


A  2  1 1 


1 1
3 
3
 
b  5
 
8
 x1  2 x2  4 x3  3

 2 x1  x2  x3  5

 x1  x2  3x3  8
 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
1 2  4


A1  0  5 9 


0  1 7 
3
 
b1   1
 
 5 
 x1  2 x2  4 x3  3

  5 x2  9 x3  1

  x2  7 x3  5
Esempio (2)
 Seconda iterazione:
 facciamo in modo che a32=0
 terza riga  terza riga – (a32/a22)seconda riga
1 2  4


A1  0  5 9 


0  1 7 
4 
1 2


A2  0  5
9 


0 0 26 / 5
3
 
b1   1
 
 5 
 3 


b2   1 


26 / 5
 x1  2 x2  4 x3  3

  5 x2  9 x3  1

  x2  7 x3  5
x  2x  4x  3
2
3
 1

  5 x2  9 x3  1
 26
26
x


3
5
5

Esempio (3)
 Il sistema può quindi essere risolto per sostituzione
partendo dall’ultima equazione:
x3  1
9 x3  1
x2 
2
5
x1  3  2 x2  4 x3  3
Considerazioni sul metodo di Gauss
 Anche 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-Jordan
 negli 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 inferiore
 U è una matrice triangolare superiore
 a11

 a21

A  LU  



a
 N1
a12
a13
a22
a23

aN 2
aN 3
 a1N   11
0
0  0   11 12
 
 
 a2 N   21  22
0  0   0  22
 
 




 
 
 
 
 
 
 a NN   N 1  N 2  N 3   NN   0
0
13  1N 
 23
0

 2N 





  NN 
 Se tale decomposizione è possibile, il sistema Ax=b si
trasforma nel modo seguente:
Ax  b  LU x  b  LUx  b
Decomposizione LU (2)
 Per risolvere il sistema di partenza occorre dunque risolvere
i due sistemi:
 Ly = b
 Ux = y
 Il vantaggio di questo approccio è che entrambi i sistemi
hanno una matrice dei coefficienti di forma triangolare, e
possono essere risolti per sostituzione:
y1 
b1
11
i 1

1 
yi 
b


y
 i  ij j 
 ii 
j 1

xN 
i  2,3,, N
yN
 NN
N

1 
xi 
 yi    ij x j 
 ii 
j i 1

i  N  1, N  2,,1
Decomposizione LU (3)
 Se vale la decomposizione A=LU si ha, dalla regola del
prodotto:
N
aij    ik  kj  ai11 j ai 2  2 j    aiN  Nj
k 1
 in totale si possono scrivere N2 equazioni di questo
tipo
 ciascuna 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 incognite
 per effettuare la decomposizione occorre fissare N
incognite
Algoritmo di Crout (1)
 Si fissano i termini diagonali della matrice L ponendo:
 ii  1
i  1,2,, N
 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:
aij   i11 j   i 2  2 j     ii  ij

N

aij    ik  kj  aij   i11 j   i 2  2 j     ii  jj
k 1

aij   i11 j   i 2  2 j     ij  jj
i j
i j
i j
 Nell’algoritmo di Crout si sfruttano queste relazioni
per ricavare, in maniera iterativa, i coefficienti ij e ij
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:
i 1
aij   i11 j   i 2  2 j     ii  ij    ik  kj   ij 
k 1
i 1
 ij  aij    ik  kj
se i  1 la sommatoria
è da considerar si nulla

k 1
 Facciamo ora variare l’indice di riga i nell’intervallo
i=j+1,j+2,...,N:
j 1
aij   i11 j   i 2  2 j     ij  jj    ik  kj   ij  jj 
k 1
j 1

1 
 ij 
a



(se j  1 la sommatoria è da considerar si nulla)

ij
ik kj 

 jj 
k 1

 I coefficienti che compaiono nelle equazioni sono stati
determinati nelle iterazioni precedenti
Esempio (1)
 Applichiamo l’algoritmo di Crout per decomporre la
matrice:
1 2  4 


A  2  1 1 


1 1
3 
 Alla fine del procedimento determineremo due matrici
L e U tali che:
A  LU
0 0
1


L   21 1 0


 31  32 1
 11 12

U   0  22

 0
0
13 

 23 

 33 
Esempio (2)
 L’algoritmo di Crout viene applicato colonna per
colonna
 si determinano prima i valori di ij e poi quelli di ij
 Colonna 1 (j=1):
11  a11  1
 21 
 31 
1
11
1
11
a21   1  2  2
1
a31   1 1  1
1
 Colonna 2 (j=2):
12  a12  2
 22  a22   2112  1  2  2  5
 32 
1
 22
a32   3112  
1
1
 1  1 2 
5
5
Esempio (3)
 Colonna 3 (j=3):
13  a13  4
 23  a23   2113  1  2   4  9
 33  a33   3113   32  23  3  3 4   9 
1
5
26
5
 Riassumendo, si ha:
1 0 0 


L   2 1 0


1 1 / 5 1
4 
1 2


U  0  5
9 


0 0 26 / 5
e si può verificare che A=LU
Commenti sull’algoritmo di Crout
 Anche 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 LU
 Calcolo 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
 N
 N
det A  det L  det U     ii    ii     ii
 i 1  i 1  i 1
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:
 1  1 0

2  2  2

 0 3 3
A





0 0 0


0

0

0
 N
0

0

0





 N 
 In generale gli elementi sulla diagonale principale sono
più grandi rispetto a quelli fuori diagonale:
i  i   i
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 normale
 non è 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 sostituzione
 poiché 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:
5 x1  x2  7
5  1 0 0 



 x1  5 x2  2 x3  15
1 5  2 0 

A

0  1  8 3 
 x2  8 x3  3x4  8



0 0
3x3  10 x4  13
3 10
7
 
 15 
B 
 8
 
 13 
la cui soluzione è x1=2, x2=3, x3=1, x4=1
 Risolviamo il sistema trasformando la matrice A in una
matrice triangolare con il metodo di eliminazione di
Gauss
Esempio (2)
 Prima iterazione:
 seconda riga  seconda riga - 2/1prima riga
5  1 0 0 


1 5  2 0 

A
0  1  8 3 


0 0
3 10
0 0
5  1


0 26 / 5  2 0 

A1  
0  1  8 3 


0
0
3 10
7
 
 15 
B 
 8
 
 13 
 7 


68 / 5

B1  
 8 


 13 
Esempio (3)
 Seconda iterazione:
 terza riga  terza riga - 3/2seconda riga
0 0
5  1


0 26 / 5  2 0 

A1  
0  1  8 3 


0
0
3 10
0
0
5  1


0 26 / 5
2
0

A2  
0
0
 109 / 13 3 


0
0
3
10
 7 


68 / 5

B1  
 8 


 13 
 7 


 68 / 5 

B2  
 70 / 13


 13 
Esempio (4)
 Terza iterazione:
 quarta riga  quarta riga - 4/3terza riga
0
0
5  1


0 26 / 5
2
0

A2  
0
0
 109 / 13 3 


0
0
3
10
 7 


 68 / 5 

B2  
 70 / 13


 13 
0
0
5  1



0 26 / 5

2
0

A3  
0

0
 109 / 13
3


0
0
0
1207 / 109
7




 68 / 5 

B3  
  70 / 13 


1207 / 109
Esempio (5)
 A questo punto il sistema può essere risolto per
sostituzione:
1207
1207
x4 
 x4  1
109
109
109
70
13  70


x3  3 x4    x3  
   3 x4   1
13
13
109  13

26
68
5  68

x2  2 x3 
 x2    2 x3   3
5
5
26  5

1
5 x1  x2  7  x1  7  x2   2
5
 In pratica, data la struttura tridiagonale, in ogni
iterazione viene cambiato soltanto un elemento della
matrice A
Metodi iterativi
 Un sistema di equazioni lineari può essere risolto in
maniera iterativa
 Si 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 sistema
 Se 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:
4x  y  z  7
4 x  8 y  z  21
 2 x  y  5 z  15
la cui soluzione è x=2, y=4, z=3
 Le equazioni del sistema possono essere riscritte nel
modo seguente:
7 yz
4
21  4 x  z
y
8
15  2 x  y
z
5
x
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:
7  yk  z k
xk 1 
4
21  4 xk  zk
yk 1 
8
15  2 xk  yk
zk 1 
5
 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
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:
 2 x  y  5 z  15
4 x  8 y  z  21
4x  y  z  7
 Ricavando x dalla prima equazione, y dalla seconda e z
dall’ultima possiamo scrivere le formule di Jacobi nel
modo seguente:
 15  y  5 z
x
2
21  4 x  z
y
8
z  7  4x  y
 15  yk  5 z k
xk 1 
2
21  4 xk  z k
yk 1 
8
z k 1  7  4 xk  yk
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
4x  y  z  7
4 x  8 y  z  21
 2 x  y  5 z  15
7 yz
x
4
21  4 x  z
y
8
15  2 x  y
z
5
 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
Metodo di Gauss-Seidel (2)
 La formula che permette di passare dalla k-esima alla
(k+1)-esima iterazione sarà dunque la seguente:
7  yk  z k
xk 1 
4
21  4 xk 1  z k
yk 1 
8
15  2 xk 1  yk 1
z k 1 
5
 Anche in questo caso si parte da una soluzione di prova
(x0,y0,z0) e si procede iterativamente
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 dominante
 Una matrice A è detta a diagonale dominante (per righe) in
senso stretto se è verificata la condizione:
N
akk   akj
k  1 N
j 1
j k
 Questo significa che in ciascuna riga l’elemento che si trova
sulla diagonale è in valore assoluto maggiore della somma
dei valori assoluti degli altri elementi
 In maniera analoga si definiscono le matrici a diagonale
dominante per colonne
 In 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
Esempio
 Consideriamo il sistema visto negli esempi precedenti:
4x  y  z  7
4 x  8 y  z  21
 2 x  y  5 z  15
 4  1 1
A   4  8 1
 2 1 5
 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:
 2 x  y  5 z  15
4 x  8 y  z  21
4x  y  z  7
 2 1 5
A   4  8 1
 4  1 1
 In questo caso, invece, la matrice A non è più strettamente a
diagonale dominante
Condizione sufficiente di convergenza
 Se 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 partenza
 Lo stesso risultato vale anche per il metodo di Gauss-Seidel
 In 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 generale
 Dato il sistema Ax=b, le formule generali utilizzate per le
iterazioni di Jacobi e di Gauss-Seidel sono le seguenti:
 Iterazione di Jacobi:
x jk 1 
b j  a j1 x1( k )    a j j 1 x (jk1)  a j j 1 x (jk1)    a jN xN( k )
a jj
N
b j   a ji xi( k )
 k 1
xj
 Iterazione di Gauss-Seidel:
x jk 1 
i 1
i j

a jj
b j  a j1 x1( k 1)    a j j 1 x (jk11)  a j j 1 x (jk1)    a jN xN( k )
a jj
j 1
x jk 1 
bj   a x
i 1
( k 1)
ji i
a jj

N
(k )
a
x
 ji i
i  j 1
Criterio di convergenza
 Per 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:
x
( k 1)
x
k 

 x
N
j 1
( k 1)
j
x

(k ) 2
j

 Condizione sulla norma generalizzata:
x
( k 1)
x
N
k 
1
  x (jk 1)  x (jk )  
j 1
 La condizione sulla norma generalizzata è più forte di
quella sulla norma euclidea