...

Metodi numerici per la risoluzione di Sistemi Lineari

by user

on
Category: Documents
39

views

Report

Comments

Transcript

Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Metodi numerici per la risoluzione di Sistemi
Lineari
Stefano Berrone
Dipartimento di Matematica
tel. 011 0907503
[email protected]
http://calvino.polito.it/~sberrone
Laboratorio di modellazione e progettazione materiali
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Problemi Triangolari - Sostituzione in avanti
Vogliamo risolvere il sistema lineare Ly = b con L triangolare
inferiore (lower triangular) non singolare.

 


b1
y1
l11 0 0 ... 0
 l21 l22 0 ... 0   y2   b2 

 


 l31 l32 l33 ... 0   y3   b3 

=


  ..   .. 
 ..
..
 .   . 
 .
.
yn
ln1 ln2 ln3 ... lnn
bn
Osservazione
Se L è non singolare, tutti gli elementi diagonali sono diversi da 0:
det(L) =
n
Y
lii
i=1
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Esplicitando le equazioni:

l11 y1
= b1 ,



 l21 y1 + l22 y2
= b2 ,
..

.



ln1 y1 + ln2 y2 + ... + lnn yn = bn .
=⇒ y1 = b1 /l11
=⇒ y2 = (b2 − l21 y1 )/l22
=⇒ yi = (bi −
i−1
X
lij yj )/lii
j=1
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodo di sostituzione in avanti (o forward substitution):
for i=1:n
y(i)=b(i);
for j=1:i-1
y(i)=y(i)-L(i,j)*y(j);
end
y(i)=y(i)/L(i,i);
end
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Un’alternativa che sfrutta la capacità di MATLAB di operare su
sottomatrici è la seguente:
y(1) =b(1)/L(1,1);
for i=2:n
y(i)=(b(i)-L(i,1:i-1)*y(1:i-1))/L(i,i);
end
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Sostituzione all’indietro
Vogliamo risolvere il sistema lineare Ux = y con U triangolare
superiore (upper triangular) non singolare.







u11 u12 u13
...
u1n
0 u22 u23
...
u2n
..
..
.
.
0
0
... un−1,n−1 un−1,n
0
0
...
0
unn







x1
x2
x3
..
.


 
 
 
=
 
 
xn
y1
y2
y3
..
.




.


yn
Osservazione
Di nuovo, dalla nonsingolarità di U segue che uii 6= 0∀i
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodi diretti
Metodi Iterativi

u11 x1 + u12 x2 +




u22 x2 +








...
...
..
.
+
+
u1n xn
u2n xn
=
=
y1 ,
y2 ,
un−1,n−1 xn−1 + un−1,n xn = yn−1 ,
unn xn
= yn .
=⇒ xn = yn /unn .
=⇒ xn−1 = (yn−1 − un−1,n xn )/un−1,n−1 .
=⇒ xi = (yi −
n
X
uij xj )/uii .
j=i+1
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodo di sostituzione all’indietro (o backward substitution):
for i=n:-1:1
x(i)=y(i);
for j=i+1:n
x(i)=x(i)-U(i,j)*x(j);
end
x(i)=x(i)/U(i,i);
end
Sfruttando le funzioni ottimizzate della libreria BLAS:
x(n)=y(n)/U(n,n);
for i=n-1:-1:1
x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);
end
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Eliminazione Gaussiana
La risoluzione di sistemi triangolari è quindi un problema
particolarmente semplice da trattare.
Idea: Nel caso di sistemi non triangolari, cercare di ricondursi ad
uno o più sistemi equivalenti di forma triangolare.
Dato un sistema lineare Ax = b le seguenti operazioni conducono
ad un sistema equivalente (ovvero con la stessa soluzione):
1
scambiare due equazioni;
2
moltiplicare un’equazione per uno scalare;
3
sostituire un’equazione con una combinazione lineare della
stessa e di un’altra equazione.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodo di eliminazione di Gauss:
consiste nell’applicare ripetutamente l’operazione 3 in modo da:
1
trasformare il sistema assegnato in un sistema a matrice
triangolare superiore
2
e quindi nel risolvere il sistema equivalente cosı̀ ottenuto con
una backward substitution.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Scriviamo il sistema Ax = b esplicitamente:

a11 x1 + a12 x2 + ... + a1n xn





 a21 x1 + a22 x2 + ... + a2n xn
a31 x1 + a32 x2 + ... + a3n xn

..


.



an1 x1 + an2 x2 + ... + ann xn
= b1 ,
= b2 ,
= b3 ,
= bn .
Supponiamo a11 6= 0
r2 ← r2 + (−
Stefano Berrone
a21
) ∗ r1
a11
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
a21 x1 +
a22 x2
+ ... +
a2n xn
=
b2
+
21
− aa11
( a11 x1 +
a12 x2
+ ... +
a1n xn
=
b1
)
0x1
(2)
(2)
(2)
(2)
+ a22 x2 + ... + a2n xn = b2 .
(2)
a22 = a22 − a12 a21 /a11 , ...,a2n = a2n − a1n a21 /a11 ,
(2)
b2 = b2 − b1 a21 /a11
Nella seconda equazione non compare più la prima incognita.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodi diretti
Metodi Iterativi
Procediamo con tutte le ultime n − 1 equazioni:
ri ← ri + mi1 ∗ r1 ,
i = 2, ..., n
ai1
mi1 = −
,
i = 2, ..., n.
a11
 (1)

a11 x1





 0
0
=⇒







0
(1)
(1)
aij = aij , bi
(2)
aij
=
(1)
aij
(1)
(1)
(1)
(2)
(2)
(2)
+ a12 x2 + ... + a1n xn = b1
(2)
(2)
(2)
+ a22 x2 + ... + a2n xn = b2
(2)
(2)
(2)
+ a32 x2 + ... + a3n xn = b3
..
.
+ an2 x2 + ... + ann xn = bn
= bi , i, j = 1, ..., n;
(1)
(2)
+ mi1 a1j , bi
(1)
= bi
Stefano Berrone
(1)
+ mi1 b1 , i, j = 2, ..., n.
Metodi numerici per la risoluzione di Sistemi Lineari
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodi diretti
Metodi Iterativi
Seconda colonna:
(2)
supponiamo a22 6= 0. Eliminiamo x2 dalle ultime n − 2 equazioni
trasformate:
(2)
ri ← ri + r2 ∗ mi2 ,
 (1)

a11 x1





 0
0
=⇒







0
(3)
(2)
mi2 = −
(1)
ai2
(2)
i = 3, ..., n
a22
(1)
(1)
(1)
(3)
(3)
(3)
+ a12 x2 + a13 x3 + ... + a1n xn = b1
(2)
(2)
(2)
(2)
+ a22 x2 + a23 x3 + ... + a2n xn = b2
(3)
(3)
(3)
+
0
+ a33 x3 + ... + a3n xn = b3
..
.
+
+ an3 x3 + ... + ann xn = bn
0
(2)
(3)
aij = aij + mi2 a2j , bi
(2)
= bi
Stefano Berrone
(2)
+ mi2 b2 , i, j = 3, ..., n
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Dopo n − 1 passi:
 (1)
(1)

a11 x1 + a12 x2



(2)

+ a22 x2

 0
0
+
0







0
+
0
(1)
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
(1)
(1)
(n)
(n)
+ a13 x3 + ... + a1n xn = b1
(2)
(2)
(2)
+ a23 x3 + ... + a2n xn = b2
(3)
(3)
(3)
+ a33 x3 + ... + a3n xn = b3
..
.
+
0
+ ... + ann xn = bn
Osservazione
Il termine noto b ha subito esattamente tutte le trasformazioni
subite dagli elementi della matrice A.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Definizione (pivot e moltiplicatori)
(i)
L’elemento a11 e i successivi elementi aii sono detti elementi
pivot. I coefficienti mij sono detti moltiplicatori.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodo di eliminazione di Gauss e fattorizzazione LU
In notazione matriciale, il k-esimo passo di eliminazione Gaussiana
puó essere rappresentato come segue:
A(k+1) = Mk A(k)
con

1

1

 mk+1,k
Mk = 

..

.
mn,k
Stefano Berrone







1
..
.
1
Metodi numerici per la risoluzione di Sistemi Lineari
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodi diretti
Metodi Iterativi
Dopo n − 1 passi:
A(n) = Mn−1 A(n−1) = Mn−1 Mn−2 A(n−2) = . . . = Mn−1 . . . M2 M1 A
|
{z
}
M
b (n) = Mn−1 . . . M2 M1 b
e
A(n) x = b (n)
con A(n) triangolare superiore
⇓
Ax = b ⇔ MAx = Mb
Posto U = MA, si può scrivere
Ux = Mb
da cui ottengo la soluzione con sostituzione all’indietro.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Esempio
Dati


4 4 8
A = 2 8 7
1 3 6
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting


12
b=9
7
calcolare la soluzione con eliminazione gaussiana.


1 0 0
M1 = − 12 1 0
− 14 0 1
e


4 4 8 | 12
M1 ∗ (A|b) = 0 6 3 | 3  .
0 2 4 | 4
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Esempio

1 0
M2 = 0 1
0 − 31
e

0
0
1


4 4 8 | 12
M2 M1 ∗ (A|b) = 0 6 3 | 3  = (U|b (3) ).
0 0 3 | 3
Con sostituzione all’indietro:
x3 = 1
x2 = (3 − 3 ∗ 1)/6 = 0
x1 = (12 − 4 ∗ 0 − 8 ∗ 1)/4 = 1
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Com’è fatta M?
- triangolare inferiore, ma

1
0
0
m21
1
0

m31 m32
1
M 6= 
 ..
..
..
 .
.
.
mn1 mn2 mn3
Poniamo

... 0
. . . 0

. . . 0

.. .. 
. .
... 1
L = M −1
U = MA =⇒ LU = LMA = A
Quindi il sistema di partenza equivale a
LUx = b.
Com’è fatta L? Suo ruolo?
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Struttura di L
L triangolare inferiore.
Il calcolo è meno complicato e costoso di quanto sembri a prima
vista: non occorre costruire M e poi calcolarne l’inversa.


1
0
0
...
0
−m21
1
0
...
0



...
0
L = −m31 −m32 1



..
.
.


.
.
−mn1 −mn2 . . . −mn,n−1 1
Osservazione
L può quindi essere costruita durante l’eliminazione gaussiana
semplicemente conservando gli elementi mij e cambiandogli il
segno.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodi diretti
Metodi Iterativi
Fattorizzazione LU
Abbiamo fattorizzato A = LU con
L triangolare inferiore
U triangolare superiore
Utilizzo per risolvere il sistema lineare:
Ax = b
LUx = b
L |{z}
Ux = b
y
⇓
Ly = b,
Stefano Berrone
Ux = y
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Il problema di partenza può essere affrontato in due fasi, risolvendo
due sistemi a matrice triangolare:
1
prima determiniamo y risolvendo Ly = b;
2
una volta noto y , determiniamo x risolvendo Ux = y .
Osservazione
La risoluzione del sistema lineare Ly = b è del tutto equivalente ad
effettuare sul vettore b le trasformazioni dell’eliminazione
gaussiana. Quindi il vettore y calcolato come soluzione del sistema
lineare Ly = b coincide con b (n) !
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Meglio fattorizzazione LU + 2 sistemi lineari triangolari, o
eliminazione gaussiana senza costruire esplicitamente L?
Con l’eliminazione Gaussiana, b subisce le stesse trasformazioni di
A. Su un solo sistema sono equivalenti. Se devo risolvere un
secondo sistema lineare
Ax = c
con eliminazione gaussiana devo rieseguire tutto il
procedimento per effettuare le trasformazioni su c e ottenere
c (n) ;
con fattorizzazione LU già calcolata devo solo risolvere i
sistemi Ly = c e Ux = y .
È veramente più conveniente fare cosı̀ piuttosto che rieffettuare
ex-novo l’eliminazione gaussiana sul sistema Ax = c?
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Costo computazionale
Per costo computazionale di un algoritmo intendiamo l’ordine di
grandezza del numero di moltiplicazioni o divisioni effettuate,
calcolato in funzione della dimensione n del problema (in questo
caso l’ordine della matrice).
Costo computazionale per calcolo fattorizzazione LU
(coincidente con quello della sola eliminazione gaussiana):
3
n
O
moltiplicazioni.
3
Il costo computazionale soluzione di un sistema triangolare
(superiore o inferiore indifferentemente):
2
n
O
moltiplicazioni.
2
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Risolvendo un sistema lineare con fattorizzazione LU la maggior
parte del costo risiede nella fattorizzazione e non nella soluzione
dei sistemi triangolari ottenuti: il costo totale è infatti di
O(
n3
n3
) + O(n2 ) ' O( )
3
3
operazioni, essendo (per n grande) n2 trascurabile rispetto a n3 .
Osservazione
I metodi di fattorizzazione sono quindi particolarmente convenienti
quando si debbano risolvere un certo numero di sistemi lineari in
cui la matrice dei coefficienti è sempre la stessa e cambia solo il
vettore dei termini noti.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Esempio
Data la matrice
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting


4 4 8
A = 2 8 7 
1 3 6
calcolarne la fattorizzazione LU.
Basta recuperare i moltiplicatori:


4 4 8
U = 0 6 3
0 0 2

1
L =  12
1
4
Osservazione: det(A) = det(L) · det(U) = 1 ·
Stefano Berrone

0 0
1 0
1
3 1
Qn
i=1 uii
= 48.
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Osservazione
L può essere memorizzata nella stessa matrice A via via che in
questa vengono azzerati gli elementi sotto la diagonale.





4 4 8
4
4 4 8
2 8 7 −→  1 6 3 −→  1
2
2
1
1
1 3 6
4 2 4
4
Stefano Berrone

4 8
6 3
1
3 2
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
L’eliminazione di Gauss / fattorizzazione LU può essere sempre
portata a termine con successo?
Definizione
Sia A ∈ Rn×n ; definiamo Ak , con k = 1, ..., n, sottomatrice
principale di testa di ordine k la matrice ottenuta intersecando le
prime k righe e k colonne di A.
Teorema
Sia A una matrice di ordine n. Se Ak è non singolare per
k = 1, ..., n − 1 allora esiste ed è unica la fattorizzazione LU di A.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Osservazione
Le matrici a diagonale dominante per righe o per colonne e le
matrici simmetriche definite positive soddisfano le ipotesi del
teorema precedente.
Osservazione
Il teorema precedente non esclude che la matrice A sia singolare.
Se anche A è singolare, ma det(Ak ) 6= 0, k = 1, ..., n − 1, esiste
comunque la fattorizzazione LU, con la matrice U che avrà un
elemento uii nullo.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
La Strategia del Pivoting
Esempio
Che sucede se voglio applicare eliminazione gaussiana a
0 1
A=
?
1 0
Occorre ricorrere alle matrici di permutazione.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Definizione
Sia Pij la matrice ottenuta dalla matrice identità permutando la
i-esima e la j-esima riga:

1 0 0
0 1 0

0 0 1




i
0 0 0
Pij = 



j 
0 0 0


0 0 0
...
...
...
...
...
...
i
0
0
0
..
.
0
..
.
1
..
.
0
...
...
...
...
...
...
j
0
0
0
..
.
1
..
.
0
..
.
0
...
...
...
...
...
...

0
0

0




0 .



0



1
Le matrici Pij sono dette matrici di permutazione semplice.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Data una matrice A
la moltiplicazione a sinistra di A per Pij (Pij A) produce lo
scambio della i-esima e j-esima riga di A;
la moltiplicazione a destra di A per Pij (APij ) produce lo
scambio della i-esima e j-esima colonna di A.
Se moltiplichiamo più matrici di permutazione semplice fra loro,
otteniamo una matrice (che indicheremo con P e che chiameremo
matrice di permutazione) il cui effetto su una matrice A sarà
quello di effettuare tutti gli scambi di righe (colonne) associate alle
matrici di permutazione semplice.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodi diretti
Metodi Iterativi
Teorema
Sia A una matrice di ordine n. Allora esiste una matrice di
permutazione P per cui si può ottenere la fattorizzazione LU di A,
cioè:
PA = LU.
Si osservi che, data A, vi possono essere diverse matrici di
permutazione P tali che PA ammette fattorizzazione LU.
Come individuare una matrice di permutazione?
⇓
Strategia di pivoting parziale.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Pivoting Parziale
Se durante l’eliminazione di Gauss al passo k-esimo risulta
(k)
akk = 0, scambiamo due equazioni in modo che il nuovo
(k)
elemento pivot akk = 0 sia diverso da zero.
La scelta è sempre tra le equazioni che seguono la k-esima.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Esempio
Consideriamo il seguente sistema lineare:

 x1 + x2 + x3 = 1
x1 + x2 + 2x3 = 2

x1 + 2x2 + 2x3 = 1
Eliminazione gaussiana: m21 = −1 e m31 = −1 =⇒

 x1 + x2 + x3 = 1
0 + 0 + x3 = 1

0 + x2 + x3 = 0
(2)
a22 = 0!
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Esempio (segue)

 x1 + x2 + x3 = 1
0 + x2 + x3 = 0
;

0 + 0 + x3 = 1
Questo è già in forma triangolare e si può ottenere per sostituzione
all’indietro la soluzione:

 x3 = 1
x2 = −x3 = −1

x1 = 1 − x2 − x3 = 1
Abbiamo det(A) = − det(U) = −1; in generale,
det(A) = (−1)s · det(U), dove s è il numero di permutazioni di
riga effettuate.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Fattorizzazione PA = LU
Cosa accade alla fattorizzazione LU con il pivoting?
A(k+1) = Mk Pk A(k)
n − 1 passi:
A(n) = Mn−1 Pn−1 A(n−1) = ... = Mn−1 Pn−1 Mn−2 Pn−2 . . . M2 P2 M1 P1 A
U = Mn−1 Pn−1 Mn−2 Pn−2 Mn−3 . . . M2 P2 M1 P1 A
= Mn−1 M̄n−2 Pn−1 Pn−2 Mn−3 . . . M2 P2 M1 P1 A
= Mn−1 M̄n−2 . . . M̄2 M̄1 Pn−1 Pn−2 . . . P2 P1 A
|
{z
}|
{z
}
M
P
= MPA
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
U = MPA
M −1 U = PA
Posto L = M −1 abbiamo determinato fattorizzazione LU di PA,
i.e. fattorizzazione PA = LU.
Fattorizzazione PA = LU per risolvere un sistema lineare:
1
Ly = Pb
2
Ux = y
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
In generale la strategia di permutare equazioni/righe viene usata
anche se si incontrano elementi pivot diversi da zero per
aumentare la stabilità dell’eliminazione gaussiana.
Si può infatti dimostrare che per aumentare la stabilità è opportuno
evitare moltiplicatori troppo grandi, poiché amplificano gli errori
di arrotondamento, con conseguente instabilità dell’algoritmo.
Per prevenire la crescita dei moltiplicatori è necessario prevenire
elementi pivot troppo piccoli.
Per raggiungere questo obiettivo possiamo applicare la strategia
del pivoting parziale ammettendo scambi di righe anche se al
(k)
passo k l’elemento pivot akk non è nullo.
Più precisamente, al passo k si determina l’elemento di modulo
(k)
massimo fra tutti gli elementi ark , con r ≥ k; una volta
individuato l’indice r corrispondente all’elemento di modulo
massimo si scambiano fra loro la riga k-esima e r -esima.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
In questo modo ci si assicura che ad ogni passo si sta usando
l’elemento pivot più grande possibile. Questa strategia
incrementa in modo significativo la stabilità dell’algoritmo.
Osservazione
Se la matrice A è a diagonale dominante per colonne è possibile
dimostrare che la strategia del pivoting parziale è del tutto
superflua in quanto non provoca alcun scambio di
equazioni.Inoltre, se la matrice A è simmetrica e definita positiva, il
metodo di Gauss senza pivoting risulta numericamente stabile,
quindi non c’è bisogno di effettuare pivoting. Per tale classe di
matrici esiste una diversa fattorizzazione molto più conveniente da
usare (la fattorizzazione di Cholesky).
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Come costruire L
P = Pn−1 . . . P1
U = A(n)
L =?
È facilmente ottenibile con l’eliminazione gaussiana tenendo conto
delle permutazione anche nell’immagazzinamento dei moltiplicatori
sotto la diagonale.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodi diretti
Metodi Iterativi
Esempio


1 4 8
A = 2 0 7 
4 2 6
P1 = P13 :


4 2 6
P13 A = 2 0 7
1 4 8


4 2
6
A(2) =  12 −1 4 
1
4
Stefano Berrone
7
2
13
2
Metodi numerici per la risoluzione di Sistemi Lineari
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Metodi diretti
Metodi Iterativi
Esempio (segue...)
P2 = P23 :

2
6
P23 A(2) =  14
7
2
13 
2
1
2
−1

(3)
A

1
L =  14
1
2
0
1
− 72

4
4

= 1
4
1
2


0
4 2
0 U = 0 72
1
0 0
Stefano Berrone
2
7
2
− 27
6
6
4

13 
2
41
7

13 
2
41
7
P = P23 P13


0 0 1
= 1 0 0
0 1 0
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Fattorizzazione di Cholesky
Se A è simmetrica e definita positiva, essendo l’eliminazione
gaussiana stabile senza pivoting, si può sfruttare la simmetria di A
per ridurre il costo computazionale del calcolo della fattorizzazione
usando la fattorizzazione di Cholesky anziché fattorizzazione LU.
Teorema (Fattorizzazione di Cholesky)
Sia A una matrice di ordine n simmetrica e definita positiva. Allora
esiste un’unica matrice triangolare superiore R con elementi
diagonali positivi tale che
A = R T R.
Costo computazionale per fattorizzazione di Cholesky:
O(
Stefano Berrone
n3
)
6
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Problemi triangolari
Eliminazione Gaussiana
Fattorizzazione LU
Pivoting
Comandi MATLAB
1
[L,U,P]=lu(A)
2
R=chol(A) Attenzione! chol assume che A sia simmetrica
definita positiva e usa solo la parte triangolare superiore!
3
\
Altre fattorizzazioni
1
Fattorizzazione QR (qr(A)): Q è ortogonale, R è triangolare
superiore. Molto stabile ma anche più costosa della
fattorizzazione PA = LU. Utile per altre applicazioni.
Varie
1
norm(x), norm(A): calcolo di norme
2
cond(A): numero di condizionamento
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Metodi iterativi
I metodi diretti consentono la costruzione della soluzione esatta,
a meno degli errori di arrotondamento, in un numero finito di
passi.
Per sistemi con matrici dense i metodi diretti sono generalmente
preferibili.
Per sistemi con matrici sparse e di ordine elevato, i metodi diretti
diventano impraticabili a causa del fill in ⇒ non può essere
sfruttata la sparsità della matrice.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
I metodi iterativi lasciano inalterata la struttura di A e
costruiscono una successione infinita di vettori {x (k) }k≥0 , che
sotto opportune ipotesi converge alla soluzione x.
Per i metodi iterativi dovremo discutere:
1
condizioni di applicabilità
2
costo computazionale (di ogni iterazione)
3
problemi connessi alla convergenza: convergenza si/no e
velocità di convergenza.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Definizione
Una successione di vettori {x (k) }k≥0 di Rn si dice convergente ad
un vettore x ∈ Rn se esiste una norma per cui
lim kx (k) − xk = 0,
k→∞
ed in tal caso si pone
lim x (k) = x.
k→∞
Osservazione
Grazie all’equivalenza delle norme su Rn tale definizione risulta
indipendente dalla norma scelta in Rn ; se una successione di vettori
converge in una norma di Rn , allora converge in tutte le norme di
Rn .
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Osservazione
La condizione di convergenza in una qualunque norma di Rn si
traduce in una condizione di convergenza per componenti:
lim x (k) = x
k→∞
⇐⇒
(k)
lim x
k→∞ i
= xi , ∀i = 1, ..., n.
Simili argomenti valgono anche per successioni di matrici
{A(k) }k≥0 . Inoltre vale il seguente teorema:
Teorema
Sia A ∈ Rn×n . Allora indicando con ρ(A) il raggio spettrale di A
lim Ak = 0 ∈ Rn×n ⇐⇒ ρ(A) < 1
k→∞
Ricordiamo che per ogni norma di matrice compatibile con una
norma di vettore si ha ρ(A) ≤ kAk.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Generalità sui metodi iterativi
Idea base: data una stima iniziale x (0) della soluzione del sistema
lineare Ax = b, con A ∈ Rn×n , det(A) 6= 0, b ∈ Rn , si costruisce
una successione di vettori {x (k) }k≥0 che converga (si spera) ad x,
risolvendo ad ogni iterazione dei sistemi lineari molto più semplici
di quello di partenza.
Come?
Consideriamo uno splitting della matrice A del tipo
A = M + N,
det(M) 6= 0.
Allora si ha
Ax = (M + N) x = b ⇐⇒ Mx = −Nx + b.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Supponiamo di essere all’iterazione k + 1, di avere quindi a
disposizione x (k) e di voler calcolare x (k+1) .
Definiamo x (k+1) in modo che valga
Mx (k+1) = −Nx (k) + b
Allora abbiamo
x (k+1) = −M −1 Nx (k) + M −1 b
che può essere scritto in forma compatta come
x (k+1) = Bx (k) + c,
avendo posto B = −M −1 N e c = M −1 b. La matrice
B = −M −1 N
è la matrice d’iterazione ed individua il particolare metodo; è
fondamentale nello stabilire le condizioni e la velocità di
convergenza.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Proprietà (consistenza)
Se x (k) → x̃, allora x̃ soddisfa l’equazione
x̃ = B x̃ + c
e quindi il limite della successione x̃ è proprio la soluzione x di
Mx = −Nx + b (cioè di Ax = b).
Osservazione
Al variare della stima iniziale x (0) si ottengono diverse successioni
{x (k) }k≥0 .
Definizione
Un metodo iterativo si dice convergente se per ogni stima iniziale
x (0) la successione {x (k) }k≥0 è convergente.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Studio della convergenza dei metodi iterativi
Sotto quali condizioni la successione {x (k) }k≥0 generata da un
metodo iterativo converge?
Se la successione converge ad un vettore x, questo è senz’altro
soluzione del sistema lineare. Definiamo l’errore e (k) commesso al
passo k come la differenza fra l’approssimazione corrente e la
soluzione:
e (k) = x (k) − x.
Chiaramente, x (k) → x ⇐⇒ e (k) → 0.
Studiamo l’errore:
e (k) = x (k) − x = Bx (k−1) + c − x
= Bx (k−1) + c − (Bx + c) = B(x (k−1) − x) = Be (k−1) .
Allora
e (k) = Be (k−1) = B(Be (k−2) ) = B 2 e (k−2) = ... = B k e (1) = B k e (0) .
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Poiché
e (k) = B k e (0)
abbiamo
lim x (k) = x ⇔ lim e (k) = 0 ⇔ lim B k e (0) = 0 ⇔ lim B k = 0
k→∞
k→∞
k→∞
k→∞
Quindi c’è convergenza se e solo se ρ(B) < 1.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Teorema
Il metodo iterativo definito dalla matrice di iterazione B converge
per ogni x (0) ∈ Rn se e solo se
ρ(B) < 1.
La condizione necessaria e sufficiente espressa dal teorema per la
convergenza non è di facile verifica. Una condizione sufficiente è
espressa dal seguente Teorema.
Teorema
Sia k · k una norma di matrice naturale. Se kBk < 1 il metodo
iterativo converge.
Dimostrazione: È sufficiente ricordare che per ogni norma naturale
ρ(B) ≤ kBk
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Il raggio spettrale indica non solo se il metodo iterativo
converge ma anche quanto velocemente converge.
Intuitivamente, la velocitá di convergenza è una misura di
quanto velocemente l’errore tende a zero.
La velocitá di convergenza di un metodo iterativo è tanto
maggiore quanto più il raggio spettrale della matrice di
iterazione è vicino a 0; viceversa è tanto più lenta quanto più
è vicino a 1.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Metodo di Jacobi
Detto anche metodo delle sostituzioni simultanee.
Consideriamo il seguente splitting di A:


F

D
A = E + D + F,
A=
E
dove
E è la parte di A sotto la diagonale,
2 D è la diagonale di A,
3 F è la parte di A sopra la diagonale
Prendiamo M = D, N = E + F . Matrice d’iterazione BJ :
1
BJ = −D −1 (E + F )
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Data l’iterata corrente x (k) ∈ Rn , la x (k+1) ∈ Rn viene calcolata
tramite il seguente algoritmo:
for i = 1 : n P
Pn
(k+1)
(k)
(k)
xi
= bi − i−1
a
x
−
a
x
/aii
j=1 ij j
j=i+1 ij j
end
In forma compatta, può essere scritto
x (k+1) = D −1 b − Ex (k) − Fx (k) .
Osservazione
Se aii = 0 per qualche i, D non è invertibile. Si possono però
preliminarmente riordinare le equazioni del sistema lineare in modo
da portare in posizione ii un elemento non nullo. Ciò è sempre
possibile se A è non singolare.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Metodo di Gauss-Seidel
Consideriamo sempre lo splitting
A = E + D + F,
ma M = E + D e N = F .
La matrice di iterazione:
BGS = −(E + D)−1 F
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Data l’iterata corrente x (k) ∈ Rn , la x (k+1) ∈ Rn viene calcolata
tramite il seguente algoritmo:
for i = 1 : n P
P
(k+1)
(k+1)
(k)
= bi − i−1
xi
− nj=i+1 aij xj
/aii
j=1 aij xj
end
In forma compatta, può essere scritto
x (k+1) = D −1 b − Ex (k+1) − Fx (k) .
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Test di arresto
Un metodo iterativo fornisce una successione infinita di iterate x (k) .
Questo diventerà un algoritmo solo se stabiliremo in modo preciso
quando interrompere il procedimento di generazione di iterate.
Dobbiamo ovvero stabilire un criterio di arresto delle iterazioni.
1
2
si arresta il processo iterativo quando risulta
(k+1)
x
− x (k) ≤ tol,
x (k+1) dove tol è una tolleranza relativa fissata sulla base di
opportuni criteri di utilizzo della soluzione.
oppure se
(k+1)
(k) − x ≤ tol,
x
dove tol è questa volta una tolleranza assoluta.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Una strada alternativa per definire un test d’arresto consiste nel
controllo del residuo dell’equazione. Sia x (k+1) la soluzione
ottenuta al passo k + 1, il residuo corrispondente è definito come
r (k+1) = b − Ax (k+1) .
Test di arresto relativo sul residuo:
b − Ax (k+1) ≤ tol.
kbk
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Osservazione
Residuo ”piccolo” non garantisce errore ”piccolo”.
(k) (k)
(k)
−1
−1 (k)
e
=
x
−
x
=
x
−
A
b
=
A
Ax
−
b
= −A−1 r (k) ≤ A−1 r (k) .
Ricordando
1
kx k
≤
kAk
kbk
(v. conti su condizionamento) abbiamo
(k) e (k) r −1 (k) k A k
≤ A
= K (A)
r kx k
kbk
kbk
Quindi
(k) r kbk
≤ tol
=⇒
(k) e kx k
Stefano Berrone
≤ K (A)
(k) r kbk
≤ K (A)tol
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Teorema (Condizioni di convergenza per il metodo di Jacobi)
Per matrici a predominanza diagonale stretta (per righe o per
colonne indifferentemente) il metodo di Jacobi converge.
Teorema (Condizioni di convergenza per il metodo di Gauss-Seidel)
Per matrici a predominanza diagonale stretta (per righe o per
colonne) il metodo di Gauss-Seidel converge.
Teorema (Condizioni di convergenza per il metodo di Gauss-Seidel)
Per matrici simmetriche definite positive il metodo di Gauss-Seidel
converge.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Confronto delle proprietà di convergenza del metodo di
Jacobi e di Gauss-Sidel
In generale la convergenza di uno dei metodi non implica la
convergenza dell’altro e viceversa.
Tuttavia quando entrambi convergono la velocità di convergenza
del metodo di Gauss-Seidel è generalmente superiore.
Teorema (Teorema di Stein-Rosenberg)
Sia A ∈ Rn×n con aij ≤ 0, ∀i 6= j e aii > 0, allora si verifica uno e
uno solo dei seguenti risultati:
1
0 < ρ(BGS ) < ρ(BJ ) < 1,
2
1 < ρ(BJ ) < ρ(BGS ),
3
ρ(BGS ) = ρ(BJ ) = 0,
4
ρ(BGS ) = ρ(BJ ) = 1.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Metodi diretti
Metodi Iterativi
Generalità
Metodo di Jacobi
Metodo di Gauss-Seidel
Test di arresto
Condizioni di convergenza per Jacobi e Gauss-Seidel
Teorema
Sia A ∈ Rn×n una matrice tridiagonale con elementi diagonali non
nulli, allora
ρ(BGS ) = ρ2 (BJ ),
cioè il metodo di Gauss-Seidel ed il metodo di Jacobi convergono o
divergono simultaneamente e il tasso asintotico di convergenza del
metodo di Gauss-Seidel è doppio di quello del metodo di Jacobi.
Stefano Berrone
Metodi numerici per la risoluzione di Sistemi Lineari
Fly UP