Risoluzione di Sistemi lineari - dm.unibo.itguerrini/html/an_08_09/sistemi_lineari_completi.pdf ·...

Carla Guerrini 2008/12/09

Risoluzione di Sistemi lineari

La risoluzione di un sistema lineare con matrice dei coefficienti uguale a una matricediagonale produce sicuramente l’algoritmo piu semplice. Sia d = [d1, d2 . . . , dn] ladiagonale principale e b = [b1, b2, . . . , bn] il termine noto, se tutti i di sono diversi da zero

la soluzione del sistema corrispondente esiste ed e unica e xi =dibi

se un termine di e nullo

e il corrispondente bi e diverso da zero allora non ci sono soluzioni, mentre se un di = 0 e ilcorrispondente bi = 0 si hanno infinite soluzioni. La figura sottostante mostra rispettivamenteuna matrice triangolare superiore, una matrice triangolare iunferiore e una matrice in cui sonoevidenziate le diagonali.

diagonaliaij=0

per i>j

aij=0

per i<j

Si consideri per esempio un sistema 3× 3 triangolare superiore

0@

u11 u12 u130 u22 u230 0 u33

1A0@

x1x2x3

1A =

0@

b1b2b3

1A

1


Se gli elementi sulla diagonale principale sono tutti diversi da zero, allora la soluzione puoessere calcolata partendo da x3

x3 =b3

u33

x2 =b2−u23x3

u22

x1 =b1−u12x2−u13x3

u11

Se il sistema ha dimensione n la riga i-esima sara

uiixi + uii+1xi+1 + . . . + uinxn = bi

da cui

xi =bi −

Pnk=i+1 uikxk

uiiquesto calcolo deve essere effettuato per i = n : −1 : 1. Quindi una possibileimplementazione MatLab per un sistema di dimensione n Ux = b puo essere

%in forma vettorialex(n)=b(n)/U(n,n); x(n)=b(n)/U(n,n);for i=n-1:-1:1 for i=n-1:-1:1

for j=i+1:n x(i)=b(i)-U(i,i+1:n)*x(i+1:n);b(i)=b(i)-U(i,j)*x(j); x(i)=x(i)/U(i,i);

end endx(i)=b(i)/U(i,i);

end

La struttura dell’algoritmo e triangolare. Per generare x(i) si compiono (i-1) moltiplicazioni,1 divisione e i-1 sottrazioni complessivamente 2i-1 operazioni. Il costo computazionale dellasostituzione all’indietro quindi e dato da

1 +nX

i=2

(2i− 1) = n2

2


Se La matrice dei coefficienti e triangolare inferiore

L =

0BBBBBB@

`1,1`2,1 `2,2`3,1 `3,2 1...

. . . . . .

`n,1 . . . ln,n−1 `n,n

1CCCCCCA

0BBB@

x1x2

x3...

xn

1CCCA =

0BBBBB@

b1b2b3...bn

1CCCCCA

la soluzione si ottiene con la sostituzione in avanti, mediante il seguente algoritmo:

%in forma vettorialex(1)=b(1)/L(1,1); x(1)=b(1)/U(1,1);for i=2:n for i=2:n

for j=1:i-1 x(i)=b(i)-L(i,1:i-1)*x(1:i-1);b(i)=b(i)-L(i,j)*x(j); x(i)=x(i)/L(i,i);

end endx(i)=b(i)/L(i,i);

end

Anche qui la struttura dell’algoritmo e triangolare. Per generare x(i) si compiono 2i-1 operazioni, rispettivam,ente i-1 moltiplicazioni, 1 divisione e i-1 sottrazioni. Il costocomputazionale della sostituzione in avanti e dato da

1 +nX

i=2

(2i− 1) = n2

Se L e una matrice triangolare bassa con elementi diagonali unitari, allora l’algoritmo non hadivisioni,

x(1)=b(1)for i=2:nb(i)=b(i)-L(i,1:i-1)*b(1:i-1)

end

3


Esempio

0@

`11 0 0`21 `22 0`31 `32 `33

1A0@

x1x2x3

1A =

0@

b1b2b3

1A

x1 =b1`11

x2 =b2−`21x1

`22

x3 =b1−`31x1−`32x2

`33

1 operazione3 operazioni

5 operzioni

In conclusione se troviamo la fattorizzazione LU della matrice A dove L e triangolare inferioree U e triangolare superiore, cioe

An×n = Ln×nUn×n

allora possiamo facilmente risolvere il sistema osservando che

Ax = b ⇒ LUx = b

ponendo Ux = y,si puo risolvere il sistema

Ly = b

per sostituzione in avanti, una volta determinato y si risolve per sostituzione all’indietro ilsistema

Ux = y.

4


Sistemi Lineari

Nelle scienze applicate la soluzione di problemi anche complessi viene spesso ricondottaalla risoluzione di uno o piu sistemi lineari della forma

Ax = b

dove A e una matrice quadrata di dimensione n × n di elementi ai,j , reali o complessi,mentre x e b sono vettori colonna di dimensione n che rappresentano rispettivamente il vettoresoluzione incognito e il vettore termine noto. Il sistema puo essere riscritto per compeneticome segue: 8

>>><>>>:

a1,1x1 + a1,2x2 + . . . + a1,nxn = b1a2,1x1 + a2,2x2 + . . . + a2,nxn = b2......

...an,1x1 + an,2x2 + . . . + an,nxn = bn

Per esempio

Si vuol determinare la situazione di equilibrio tra la domanda e l’offerta di determinatibeni. In particolare consideriamo n industrie che producono n prodotti diversi. Esse dovrannofar fronte alla domanda interna (ossia alle quantita di beni necessarie alle industrie stesse perprodurre) ed alla domada esterna dovuta ai consumatori. Indichiamo con xi, i = 1 . . . , nnumero di unita del bene i−esimo prodotte dall’industria i−esima e con bi, i = 1, . . . , nil numero di unita del bene i-esimo assorbite dal mercato. Infine denotiamo con ci,j la frazionedi xi necessaria all’industria j − esima per la produzione del prodotto j − esimo. Sesupponiamo che le funzioni di trasformazione che legano i vari prodotti siano lineari ( secondoun modello proposto da Leontief), allora in una situazione di equilibrio il vettore x delleproduzioni totali dovra eguagliare la somma delle domande, ossia

x = Cx + b

5


essendo C = (ci,j) e b = (bi). Di conseguenza, in tale modello la produzione totalesoddisfa al seguente sistema lineare:

Ax = b, ove A = I − C.

b_1

b_2

b_3

c11

c12

c22

c31

c33

La soluzione del sistema esiste se e solo se la matrice A e non singolare.In linea teorica essa potrebbe essere calcolata tramite la cosidetta regola di Cramer:

xi =det(Ai)

det(A), i = 1, 2, . . . , n

6


dove Ai e la matrice che si ottiene sostituendo alla i−esima colonna di A il vettore terminenoto b e det(A) denota il determinante di A. Tale regola se si calcolano i determinanti conlo sviluppo di Laplace richiede un costo computazionale di circa 2(n + 1)!. Ad esempio su

un calcolatore in grado di eseguire 109 operazioni al secondo (cioe 1 giga-flops) sarebbero

necessarie 12 ore per risolvere un sistema di dimesione n = 15, se n = 100 circa 10143 anni.

7


Esistenza della Fattorizzazione LU di A

Sia A una matrice quadrata di ordine n, e siano Ak le sue sottomatrici principali ditesta di ordine k. Se Ak e non singolare per k = 1, 2, . . . , n− 1 allora esiste ed e unica lafattorizzazione LU di A.

A = LU,

in cui L ed U sono rispettivamente, una matrice triangolare inferiore ed una matrice triangolaresuperiore.

Le classi di matrici per cui vale il teorema sono:

• Le matrici simmetriche definite positive. Una matrice simmetrica A si dice definitapositiva se

∀x ∈ Rn

con x 6= 0 xT

Ax > 0

• le matrici a dominanza diagonale stretta per righe o per colonne. Una matrice si dice adominanza diagonale stretta per righe se

|ai,i| >X

j=1,j 6=i

|ai,j|, i = 1, 2, . . . , n,

per colonne se:

|aj,j| >X

i=1,j 6=i

|ai,j|, j = 1, 2, . . . , n,

Osserviamo che se A e non singolare anche L e U devono esserlo, cio assicura che i loroelementi diagonali siano diversi da zero.

In questo caso la risoluzione del sistema Ax = b puo essere ricondotta alla risoluzionedi due sistemi triangolari, piu precisamente:

8


Ax = b =⇒ LUx = b =⇒

Ly = b,Ux = y.

Ciascuno di essi puo essere facilmente risolto.Infatti essendo il sistema Ly = b triangolare inferiore, la prima riga del sistema sara dellaforma:

l1,1y1 = b1

da cui essendo l1,1 6= 0 si ottiene y1. Sostituendo il valore trovato per y1 nelle n-1equazioni successive si ottiene ancora un sistema triangolare inferiore, nelle sole incognitey2, y3, . . . , yn per il quale sipotra procedere in modo analogo ricavando y2. Procedendoquindi equazione per equazione calcoleremo tutte le incognite con la seguente procedura dettadi sostituzione in avanti:

8>>><>>>:

y1 = 1l1,1

b1

yi = 1li,i

“bi −

Pi−1j=1

li,jyj

”

per i = 2, . . . , n

function x=ltrisol(l,b)% function x=ltrisol(l,b)% l - matrice triangolare inferiore% b vettore dei termini noti column-vectorn=length(b);x(1,1)=b(1)/l(1,1);for i=2:n

for j=1:i-1if abs(l(i,i))<= n*epsdisp(’matrice mal condizionata’)

endb(i,1)=b(i,1)-l(i,j)*x(j,1);end

x(i,1)=b(i,1)/l(i,i);end

Vediamo quante operazioni occorrono per eseguire il codice, per calcolare le incognite

9


yi si devono effettuare i− 1 somme, i− 1 prodotti e una divisione, quindi

nX

i=1

1 + 2nX

i=1

(i− 1) = 2nX

i=1

i− n ≈ O(n2)

operazioni. Nel caso in cui si valuti la complessita computazionale in termini di solemoltiplicazioni allora

cc ≈ O(n2

2)

Con una procedura analoga potra essere risolto il sistema triangolare alto Ux = y, in questocaso la prima incognita ad essere calcolata sara xn e poi a ritroso verrannop calcolate tuttele altre incognite xi:

8>><>>:

xn = 1un,nyn

xi = 1ui,i

“yi −

Pnj=i+1 ui,jxj

”

per i = n− 1, . . . , 1

function x=utrisol(u,b)% function x=utrisol(u,b)% u - matrice triangolare superiore% b vettore dei termini noti column-vectorn=length(b);x(n,1)=b(n)/u(n,n);for i=n-1:-1:1

for j=i+1:nif abs(u(i,i))<= n*epsdisp(’matrice mal condizionata’)

endb(i,1)=b(i,1)-u(i,j)*x(j,1);end

x(i,1)=b(i,1)/u(i,i);end

10


Questo algoritmo viene chiamato della sostituzione all’indietro e richiede ancora n2

operazioni, in pratica ha la stessa complessita computazionale dell’algoritmo precedente. Quindiil problema si riduce a calcolare L ed U a partire dalla matrice A. Il procedimento costruttivoe basato sul processo di eliminazione di Gauss e permette di calcolare la fattorizzazione dellamatrice.

Metodo di Gauss Consideriamo un caso semplice per illustrare come funziona il metododi gauss per la fattorizzazione LU di una matrice A3× 3. Sia

0@

a11 a12 a13a21 a22 a23a31 a32 a33

1A

con

a11 6= 0 e det

„a11 a12a21 a22

«= a11a22 − a12a21 6= 0

Consideriamo la matrice L1

L1 =

0BB@

1 0 0

−a21a11

1 0

−a31a11

0 1

1CCA

si osservi la necessita che a11 sia diverso da zero, a11 e detto perno o pivot, allora

L1A =

0BB@

a11 a12 a13

0 a(1)22 a

(1)23

0 a(1)32 a

(1)33

1CCA

11


ove a(1)i,j

= aij −ai1a11

a1,j, i, j = 2, 3. Consideriamo ora la matrice

L2 =

0BBB@

1 0 00 1 0

0 −a(1)32

a(1)22

1

1CCCA

anche qui e necessario che a(1)22 =

a11a22−a12a21a11

6= 0; a(1)22 e detto perno o pivot.

Allora avremo

L2(L1A) =

0BB@

a11 a12 a13

0 a(1)22 a

(1)23

0 0 a(2)33

1CCA

con a(2)33 = a

(1)33 − a

(1)32

a(1)22

Chiamiamo U la matrice L2L1A si noti che L1 e L2 sono non

singolari, e quindi

A = L−11 L

−12 U

ove

L−11 =

0BB@

1 0 0a21a11

1 0a31a11

0 1

1CCA L

−12 =

0BBB@

1 0 00 1 0

0a(1)32

a(1)22

1

1CCCA

L−11 L

−12 =

0BBBB@

1 0 0a21a11

1 0

a31a11

a(1)32

a(1)22

1

1CCCCA

quest’ultima la possiamo chiamare L

EsempioSia

Ax = b con A =

0@

3 2 12 3 11 2 3

1A , e b =

0@

393426

1A

12


L1 deve essere triangolare inferiore e tale da rendere nulli gli elemnti della prima colonna diA cioe sotto l’elemento a11 = 2.

L1 =

0@

1 0 0−2/3 1 0−1/3 0 1

1A

L1A =

0@

0 0−2/3 1 0−1/3 0 1

1A0@

3 2 12 3 11 2 3

1A =

0@

3 2 10 5/3 1/30 4/3 8/3

1A

L2 deve essere triangolare inferiore e tale da rendere nulli gli elementi della seconada colonna

di L1A sotto a a(1)22 = 5/3

L2 =

0@

1 0 00 1 00 −4/5 1

1A

L2(L1A) =

0@

1 0 10 1 00 −4/5 1

1A0@

3 2 10 5/3 1/30 4/3 8/3

1A =

0@

3 2 10 5/3 1/30 0 12/5

1A = U

Allora

L2L1A = U

A = L−11 L

−12 U = LU

ove

L = L−11 L

−12 =

0@

1 0 02/3 1 01/3 4/5 1

1A

Ora avendo fattorizzato la matrice possiamo procedere alla soluzione dei sistemi triangolari

Ly = b e Ux = y

13


0@

1 0 02/3 1 01/3 4/5 1

1A0@

y1y2y3

1A =

0@

393426

1A y1 = 39;

„1 0

4/5 1

«=

„y2y3

«=

„34− 39 · 2/326− 39 · 1/3

«y2 = 8;

(1)(y3) = (13− 8 · 4/5) y3 = 33/5 = 6.6

0@

3 2 10 5/3 1/30 0 12/5

1A0@

x1x2x3

1A =

0@

398

6.6

1A x3 = 11/4 = 2.75

„3 20 5/3

«=

„x1x2

«=

„39− 11/4

8− 1/3 · 11/4

«x2 = 85/12·3/5 = 17/4 = 4.25;

(3)(x1) = 145/4− 2 ∗ 17/4

x1 = 1/3(111/4) = 111/12 = 9.25

14


Possiamo porci come obiettivo il calcolo del determinante di A quindi det(A) =det(L) ∗ det(U). Nell’esempio precedente ricaviamo che det(A) = det(L) · det(U) =1 ∗ 12 = 12

Consideriamo il caso generale con dimensione n > 3. Si osservi la seguente matrice

n× n , Li,j = I + `i,jeieTj con i > j

Li,j =

0BBBBBBBBBBB@

10 1...

. . . 1

. . . . . . . . .. . .

0... `i,j

. . . 1... . . . . . .

.... . .

0 . . . . . . . . . . . . 0 1

1CCCCCCCCCCCA

E una matrice triangolare inferiore con determinante uguale a 1. Per esempio

L2,1A =

0BBBBBB@

a1,1 a1,2 . . . a1,na2,1 + `2,1a1,1 a2,2 + `2,1a1,2 . . . a2,n + `2,1a1,na3,1 a3,2 . . . a3,n...

......

...an,1 an,2 . . . an,n

1CCCCCCA

Se a1,1 6= 0 e scegliamo

`2,1 = −a2,1

a1,1

allora verra azzerato l’elemento nella posizione (2,1), e cambieranno solamente gli elementiche appartengono alla seconda riga. Si osservi che questa moltiplicazione non ha cambiato ildeterminante della matrice di partenza. Possiamo costruire, allo stesso modo, la matrice L3,1per azzerare l’elemento di posizione (3,1), in questo caso

`3,1 = −a3,1

a1,1

15


La moltiplicazione di L3,1A portera alla costruzione di una nuova matrice in cui solo la terzariga risulta variata rispetto alla precedente. Il procedimento di annullamento procede sino afar risultare uguali a zero tutti gli elementi che appartengono alla prima colonna escluso ilprimo. Alla fine otterremo la matrice Ln,1Ln−1,1Ln−3,1 . . . L2,1A che ha lo stesso

determinente di A, la prima colonna nulla eccetto il primo elemento, che indichiamo con A(2)

i cui elementi delle altre righe sono dati da:

a(2)i,j

= ai,j + `i,1a1,j = ai,j −ai,1a1,j

a1,1

per j = 2, 3, . . . , n e i = 2, 3, . . . , n.Poniamo

L1 = Ln,1Ln−1,1Ln−3,1 . . . L2,1

si verifica facilmente che

L1 =

0BBBBB@

1`2,1 1

`3,1 0 1...

. . . . . .

`n,1 0 . . . . . . 1

1CCCCCA

L−11 =

0BBBBB@

1−`2,1 1

−`3,1 0 1...

. . . . . .

−`n,1 0 . . . . . . 1

1CCCCCA

16


La matrice A(2) = L1A ha quindi la forma

L1A =

0BBBBBBBBB@

a1,1 a1,2 . . . a1,n

0 a(2)2,2 . . . a

(2)2,n

0 a(2)3,2 . . . a

(2)3,n

......

......

0 a(2)n,2 . . . a

(2)n,n

1CCCCCCCCCA

Possiamo procedere annullando gli elementi della seconda colonna al disotto della diagonaleusando le matrici

I + `j,2ejeT2

con j > 2. Per fare questo dobbiamo supporre che a(2)2,2 6= 0. Se a

(2)2,2 = 0 possono

presentarsi due casi:

• Tutti gli elementi della seconda colonna sono nulli (ecceto il primo), allora in questo casoil determinante sarebbe nullo essendo la prima e la seconda colonna proporzionali.

• Esiste un elemento al posto k − esimo non nullo, allora si puo scambiare la seconda ela k − esima riga in modo che l’elemento nella posizione (2,2) sia non nullo.

Supponendo che a(2)2,2 6= 0, prendendo `j,2 = −

a(2)j,2

a(2)2,2

, possiamo annullare tutti gli elementi

della seconda colonna al disotto della diagonale principale. Si definisce quindi

L2 = Ln,2Ln−1,2Ln−3,2 . . . L3,2 =

0BBBBBBB@

10 10 `3,2 1

0 `4,2 0 1...

. . . . . . . . .

0 `n,2 0 . . . 1

1CCCCCCCA

17


e la nuova matrice

A(3)

=

0BBBBBBBBB@

a1,1 a1,2 a1,3 . . . a1,n

0 a(2)2,2 a

(2)2,3 . . . a

(2)2,n

0 0 a(3)3,3 . . . a

(3)3,n

......

......

...

0 0 a(3)n,3 . . . a

(3)n,n

1CCCCCCCCCA

che ha sempre lo stesso determinante di A. Si verifica facilmente che gli elementi al disottodella seconda riga sono dati da:

a(3)i,j

= a(2)i,j

+ `i,2a(2)2,j

= a(2)i,j

−a(2)i,2

a(2)2,j

a(2)2,2

Per j = 3, 4, . . . , n e i = 3, 4, . . . , n. Se il procedimento non si arresta prima,fornendoci l’informazione che il determinante e nullo si arrivano a definire le matriciLn−1, Ln−3, . . . , L2, L1 tali che:

A(n−1)

= Ln−1Ln−3 . . . L2L1A =

0BBBBBBBBB@

a1,1 a1,2 a1,3 . . . a1,n

0 a(2)2,2 a

(2)2,3 . . . a

(2)2,n

0 0 a(3)3,3 . . . a

(3)3,n

......

. . . . . ....

0 0 0 . . . 0 a(n−1)n,n

1CCCCCCCCCA

Poiche le matrici Lj sono tutte triangolari inferiori con determinante uguale a 1 la matrice

ottenuta ha lo stesso determinante di A, la nuova matrice A(n−1) e triangolare superiore

18


con determinante uguale a

a1,1a(2)2,2a

(3)3,3 . . . a

(n−1)n,n .

Se anche a(n−1)n,n risulta diverso da zero allora il determinante di A e non nullo.

Gli elementi `i,k delle matrici Li e gli elementi della matrice triangolare superiore A(n−1)

(che indicheremo poi con U) si possono ottenere con la seguente procedura Matlab:

% algoritmo di Gaussfor k=1:n-1for i=k+1:n

m(i,k)=-a(i,k)/a(k,k);for j=k+1:n

a(i,j)=a(i,j)+l(i,k)*a(k,j);endendend

Ogni matrice Li ha l’inversa che si ottiene cambiando di segno agli elementi `i,j . Seindichiamo con

L−1

= Ln−1Ln−2 . . . L2L1allora

L = (Ln−1Ln−2 . . . L2L1)−1

= L−11 L

−12 . . . L

−1n−2L

−1n−1

quindi abbiamo ottenutoA = LU

Gli elementi a(k)k,k

devono essere tutti non nulli, vengono chiamati elementi pivot. Non

e necessario memorizzare gli elementi delle matrici A(k), essi di volta in volta vengonomemorizzati al posto degli elementi di A, al termine del processo gli elementi della matrice

triangolare superiore U sono dati da ui,j = a(i)i,j

per i = 1, 2, . . . , n e j = i, i+1, . . . , n

non importa memorizzare gli elementi diagonali della matrice L, perche sappiamo che valgono1 (per costruzione).L’algoritmo appena scritto e alla base dei seguenti comandi MATLAB:

À [L, U ] = lu(A);

19


À B = inv(A);

À A\b

Costo computazionalePer ottenere la fattorizzazione LU di A occorrono le seguenti operazioni:

• Per il calcolo dei coefficienti li,j che sono in numero di n−1 nella prima colonna, n−2nella seconda ecc. fino ad 1 per la penultima, occorrono

n−1X

j=1

(n− j) =n(n− 1)

2

divisioni

• Per il calcolo degli a(k)i,k

, che sono (n− 1)2 al primo passo (k=2), (n− 2)2 al secondo

passo, fino ad 1 (k=n-1). In totale il numero degli a(k)i,k

e dato da:

n−1X

k=1

(n− k)2

=n(n− 1)(2n− 1)

6' n3

3

per ciascuno dei quali ocorrono 1 sottrazione e una moltiplicazione, quindi trascurando lepotenze piu basse, occorrono

∼ n3

3 moltiplicazioni e divisioni

∼ n3

3 sottrazioni.

Calcolo del determinante di A

La fattorizzazione A = LU ci permette di calcolare

det(A) = det(LU) = det(L)det(U) = Πni=1ui,i

20


con costo computazionale basso il tradizionale algoritmo basato sul calcolo dei determinantidei minori della matrice ha un costo computazionale dell’ordine di n!.

Calcolo del rank(A)

La fattorizzaizione A = LU ci permette di calcolare con costo computazionale basso

rank(A) = rank(U)

Calcolo dell’inversa A−1

Sia A ∈ Rn×n ivertibile, vogliamo calcolare A−1, risolvendo il sistema

AX = In

Sia

X = [x1|x2| . . . |xn],

In = [e1|e2| . . . |en],

ove In e la matrice identita di dimensione n, i vettori e1, e2, . . . , en formano la base

canonica di Rn, X ∈ Rn×n e la matrice incognita. il calcolo diretto di A−1 costa O(n4)

operazioni. Risolvere AX = In per ottenere X = A−1In e equivalente a risolvere nsistemi lineari, tutti con la stessa matrice dei coefficienti A:

8>>><>>>:

Ax1 = e1Ax2 = e2...

...Axn = en

21


Fattorizziamo A = LU una sola volta e poi risolviamo 2n sistemi lineari con matricetriangolare 8

>>>>>>>><>>>>>>>>:

(i) A = LU On3

3

(ii)

Ly1 = e1Ux1 = y1

On2

2...

...

(n)

Lyn = e1Uxn = y1

On2

2

Complessivamente il costo computazione risulta:

O(n3

3+ 2n

n2

2) = O(

4n3

3)

EsempioSupponiamo di dover risolvere il sistema

A =

24

1 1− ε 32 2 23 6 4

35

b =

24

5− ε613

35

la cui soluzione e x = (1, 1, 1)T indipendentemente dal valore di ε. Se ε assume valorimolto piccoli oppure risulta uguale a zero pur essendo A non singolare l’algoritmo scrittosopra puo comportare delle divisioni per numeri molto piccoli o addirittura per zero. L’analisidella stabilita dell’algoritmo porta alle seguenti considerazioni se l’elemento pivot al passo k,

akk,k per effetto dell’aritmetica finita risulta essere zero o un numero molto piccolo allora gli

elementi `i,k =ai,kak,k

o non sono calcolabili oppure possono diventare numeri molto grandi

e causare overflow nel calcolo degli elementi ai,k . Questo fenomeno accade quando nonvalgono le ipotesi del teorema iniziale. A tale proposito si dimostra che

22


TeoremaSia A una matrice di ordine n non singolare. Allora esiste una matrice di permutazione P percui si puo ottenere la fattorizzazione LU di PA. Cioe

PA = LU

Consideriamo le matrici di permutazione elementari Pr,s, ottenute dalla matrice identita Iscambiando la r − esima riga con la s− esima.

P =r,s

1 00 1 0

0 1

0

.

.

.

. . ..

.1001

. . 0. .

.

.

. .

s

sr

r10

0

1

11

11 ..

E facile verificare che il prodotto a sinistra di Pr,s per la matrice A e una matrice incui risultano scambiate le righe r ed s mentre le altre rimangono in variate. Nel prodotto adestra APr,s risultano scambiate le colonne.(Si osservi che det(P ) = ±1 secondo se il numero delle permutazioni e pari o dispari).Esempio

A =

24

1 2 43 5 67 8 9

35 P =

24

1 0 00 0 10 1 0

35

PA =

24

1 2 47 8 93 5 6

35 AP =

24

1 3 24 6 57 9 8

35

23


La matrice di permutazione viene determinata in modo costruttivo. Illustriamo il metodo conun esempio. Sia data la matrice

A = A0 =

24

4 −8 22 −4 61 −1 3

35

Il pivot a(0)i,1

e 4 e quindi P1 = I e

L1 =

24

1 0 0−1/2 1 0−1/4 0 1

35 ;

si costruisce

A1 = L1P1A0 =

24

4 −8 20 0 50 1 5/2

35

Il pivot a(1)i,1

e 1 le matrici P2 e L2 saranno:

P2 =

24

1 0 00 0 10 1 0

35 L2 =

24

1 0 00 1 00 0 1

35 ;

si costruisce

A2 = L2P2A1 =

24

4 −8 20 1 5/20 0 5

35 = U

Quindi abbiamo ottenutoL2P2L1P1A = U.

Si noti cheL2P2L1P1 = L2P2L1P

−12 P2P1

24


Si osservi che L2P2L1P−12 e una matrice triangolare inferiore che possiamo chiamrae T ,

se poniamo

L−1

= T , P = P2P1

otteniamo

PA = LU

Completando l’esempio abbiamo:

L = P2L−11 P

−12 L

−12 =

24

1 0 01/4 1 01/2 0 1

35

e infatti

LU =

24

1 0 01/4 1 01/2 0 1

35 ·24

4 −8 20 1 5/20 0 5

35 =

24

4 −8 21 −1 32 −4 6

35 = PA

Si noti che non si deve perder l’informazione contenuta in P nel corso dell’algoritmo, infatti sefattorizzaimo PA il sistema da risolvere sara

PAx = Pb

Metodo di Gauss con scambio delle righe e perno Massimo

Per rendere piu stabile l’algoritmo e non generare dei problemi di overflow si considerala tecnica del pivot per colonneEsempioConsideriamo la matrice

B =

»δ 11 1

–

25


ove 0 < δ < 1 e piccolo a piacere, quindi B ammette la fattorizzazione:

L =

"1 01δ

1

#U =

"δ 1

0 1− 1δ

#

In questo caso il metodo e applicabile, ma quando δ e molto piccolo per esempio δ = 10−20,supponiamo di lavorare in aritmetica finita e in doppia precisione, cioe con una unita di

macchina u =' 10−16, si ottiene

L =

"1 0

1020 1

#

U =

"10−20 1

0 1− 1020

#'"

10−20 1

0 −1020

#

poiche 1− 1020 ' −1020, se calcoliamo il prodotto

LU =

"10−20 11 0

#6="

10−20 11 1

#

Conviene percio eseguire il pivoting ad ogni passo, nel corso della costruzione della

fattorizzazione si cerca fra tutti i pivot disponibili, cioe gli elementi a(k)i,k

con i =

k, k + 1, . . . , n quello di modulo massimo.Formalmente il procedimento diventa il seguente:

• si cerca l’elemento di modulo massimo 6= 0 nella prima colonna di A supponimamo sianella r − esima riga, tale elemento deve esistere altrimenti A sarebbe singolare, siindividua cosı la matrice di permutazione P1 che scambia la riga r con la riga 1.

• si calcola la matrice L1 che azzera gli elementi sotto-diagonali ∈ alla prima colonna diP1A

26


• si cerca l’elemento di modulo massimo 6= 0 tra gli elementi della seconda colonna,escludendo la prima riga della matrice L1P1A, (se non esiste allora la matrice dipartenza o e singolare o risulta prossima a una matrice singolare), si individua cosı lamatrice di permutazione P2 ( se non si effetuano scambi di righe allora P2 = I).

• si procede ripetendo i passi precedenti...

Al termine della procedura avremo:

Ln−1Pn−1Ln−2Pn−2 . . . L2P2L1P1A = U

Si osservi che ogni matrice di scambio e involutoria cioe PiPi = P2i = I per cui si ha

immediatamente che

(P1P2 . . . Pn−2Pn−1)(Pn−1Pn−2 . . . P2P1) = I

allora l’eguaglianza precedente si puo scrivere

Ln−1Pn−1Ln−2Pn−2 . . . L2P2L1...

...(P2 . . . Pn−2Pn−1)(Pn−1Pn−2 . . . P2)P1A = U

cioe ancora posto P = Pn−1Pn−2 . . . P2P1

(Ln−1Pn−1Ln−2Pn−2 . . . L2P2L1P2 . . . Pn−2Pn−1)PA = U

Si tratta di riconoscere che la matrice tra parentesi e ancora triangolare inferiore ed invertibile,a tale proposito si osserva che la generica matrice Pi o e l’identita o scambia la riga di

indice i con una riga di indice piu elevato. E facile riconoscere che per esempio P2L1P2 eancora triangolare inferiore L1 con elementi diagonali uguali a 1, L2L1 e triangolare inferiorecon elementi diagonali uguali a 1 elementi diversi da zero solo nelle prime due colonne. Lamatrice di permutazione P3 scambia al piu la terza riga con una riga di indice piu elevato

27


e quindi si ricade sempre come sopra cioe P3L2L1P3 e dello stesso tipo, si riconoscequindi che (Ln−1Pn−1Ln−2Pn−2 . . . L2P2L1P2 . . . Pn− 2Pn− 1) e triangolareinferiore ed e invertibile in quanto i suoi termini diagonali valgono 1. Se tale matrice viene

indicata con L−1 si ottine L−1PA = U cioe

PA = LU

L’algoritmo con pivoting per colonne, diventa allora:

Per k = 1 : n− 1Per i = k + 1 : n

trova mk 3 |amk,k| = maxi=k:n |a(k)i,k|

scambio righe k ⇐⇒ mk

li,k =a(k)i,k

a(k)k,k

Per j = k + 1 : n

a(k+1)i,j

= a(k)i,j

− li,ka(k)k,j

Esempio

Consideriamo il sistema lineare Ax = b con A ∈ Rn×n di elementi

ai,j =1

i + j − 1, i, j = 1, . . . , n

Per n=5 la matrice di Hilbert assume la seguenmte forma

28


H =

26666664

1 12

13

14

15

12

13

14

15

16

13

14

15

16

17

14

15

16

17

18

15

16

17

18

19

37777775

e b scelto in modo tale che la soluzione esatta del sistema sia x = (1, 1, 1 . . . , 1)T .

La matrice introdotta e detta matrice di Hilbert ed e simmetrica definita positiva, siottiene nell’ambito delle applicazioni del metodo dei minimi quadrati. Calcoliamo in Matlabla fattorizzazione di Gauss con pivoting parziale e risolviamo i sistemi lineari al variare di n,indichiamo con xc la soluzione calcolata.

% script% test soluzione di un sistema lineare con matrice di Hilbertfor n=1:100h=hilb(n);b=sum(h’)’;xv=ones(i,1);[l,u,p]=lu(h); yc=l\(p*b);xc=u\yc;errx(n)=norm(xv-xc)/norm(xv);errR(n)=norm((l*u-p*h),inf);endfigure(1),plot(log(errx))title(’errore relativo sulla soluzioneal variare della dimensione’)figure(2),plot(errR)title(’norma della matrice residuo’)

29


0 10 20 30 40 50 60 70 80 90 100−40

−30

−20

−10

0

10

20errore relativo sulla soluzioneal variare della dimensione

0 10 20 30 40 50 60 70 80 90 1000

0.5

1

1.5

2

2.5

3

3.5

4x 10

−16 norma della matrice residuo

andamento al variare di n del logaritmo dell’errore relativo parte sinistra, a destra normainfinito della matrice residuo

Il programma Matlab lu calcola la fattorizzazione con pivoting parziale la sintassi e:

À [L, U, P ] = lu(A)

se viene richiamato nella forma [L, U ] = lu(A), non si ha la garanzia che L sia triangolareinferiore, lo e a meno dell’inversa della matrice di permutazione P che viene automaticamenteincorporata ad L.

30


Stabilita della fattorizzazione LU

Le operazioni aritmetiche che intervengono nell’algoritmo di fattorizzazione sonoeffettuate in aritmetica finita segue che questi algoritmi, invece di generare i fattori L eU per la matrice A, generano dei fattori non esatti

L = L + δL, e U = U ∗ δU

posto LU = A + δA si ha

A + δA = (L + δL)(U + δU) =

= LU + UδL + LδU + δUδL

da cui

δA = UδL + LδU + δLδU

Da cui segue che se gli elemnti di Le U sono grandi , gli elementi di δA sono grandi e quindigli errori di arrotondamento si amplificano (analisi dell’errore all’indoetro). Percio diremo chela fattorizzazione LU e stabile numericamente se gli elemnti di L e di U non sono troppograndi rispetto agli elemnti di A. In particolare se esistono delle costanti a e b indipendentida A e dall’ordine di A tali che |`ij| ≤ a e |uij| ≤ b, allora si dice che la fattorizzazioneLU e stabile in senso forte; se le costanti a e b dipendono dall’ordine di A allora si diceche la fattorizzazione LU e stabile in senso debole. Lo scopo della scelta del pivot secondoil criterio massimo e appunto quello di fare in modo che le singole compenenti delle matriciLe U siano di grandezza paragonabile a quella delle compenenti della matrice A. Osserviamoche la scelta del massimo pivot per colonne comporta che gli elementi della matrice L

|lj,k| ≤ 1 ∀j > k

⇒ ‖L‖ ≤ α ≤ n

31


quindi l’algoritmo di Gauss con pivoting parziale genera una fattorizzazione stabile in sensodebole.

EsempioData la matrice

W = (wi,j) ∈ Rn×n

,

wi,j =

8>>><>>>:

1 se i = j (ai,i = 1)

1 se i = m (ai,m = 1)

0 se i < j &i 6= m−1 se i > j

Ad esempio per m = 5

W =

266664

1 0 0 0 1−1 1 0 0 1−1 −1 1 0 1−1 −1 −1 1 1−1 −1 −1 −1 1

377775

Applicando la fattorizzazione la fattorizzazione di Gauss si ottiene:

W = LU =

266664

1 0 0 0 0−1 1 0 0 0−1 −1 1 0 0−1 −1 −1 1 0−1 −1 −1 −1 1

377775

266664

1 0 0 0 10 1 0 0 20 0 1 0 40 0 0 1 80 0 0 0 16

377775

Il fattore di crescita per la matrice W e esattamente 25−1.

Riassumendo data una matrice A ∈ Rn×n applicando la fattorizzazione di Gauss conpivoting parziale si ottiene:

‖L‖ = O(1), ‖U‖ = O(ρu)‖A‖ 6= O(u)‖A‖

33


dove il fattore di crescita ρ puo raggiungere il valore 2n−1. Questo risultato afferma chela fattorizzazione PA = LU non e sempre stabile. In questo caso parliamo di stabilita insenso debole perche dipende dal fattore di crescita ρ.

Un’altra strategia per il metodo di Gauss che consente di ridurre ulteriormente il fattore

di crescita degli elementi delle matrici A(k) e quella basata sul pivoting totale. Consistenel sostituire al k − esimo passo l’elemento nella posizione (k, k) con il massimo dellasottomatrice di ordine n − k × n − k. Occorre osservare che mentre lo scambio delle righenon affligge in alcun modo la soluzione, nel caso di scambio di colonne occorre procedere alloscambio dei corrispondenti elementi del vettore della soluzione. Con questa strategia il fattoredi crescita cresce piu lentamente che nel caso precedente ed e del tipo

g(k) =

sk ∗ Πk

j=2J

1j−1

come si vede dal grafico.

2 4 6 8 10 12 14 16 18 200

20

40

60

80

100

120

140

160

180

200f(k)=2k

g(k)=sqrt(k*prod(j(1/j−1))

34


Stabilita della Sostituzione all’indietro

Conducendo una analisi all’indietro, consideriamo la soluzione calcolata x come soluzioneesatta del problema perturabato:

(L + δL)x = y per‖δL‖‖L‖ = O(u)

Si dimostra che

‖δL‖∞ ≤ n(n + 1)

2u

analogamente per la sostituzione all’indietro.

35


Analisi del Condizionamento di unsistema Ax = b

Qualunque metodo si usi per risolvere un sistema lineare, non possiamo prescindere daglierrori di rappresentazione dei dati e dagli errori di arrotondamento nei calcoli. La soluzione diun sistema lineare, tenendo conto degli errori dovuti all’uso dell’aritmetica floating point e apossibili errori nel termine noto, puo essere considerata come soluzione esatta di un problemaperturbato del tipo:

(A + δA)(x + δx) = b + δb

ove x e la soluzione esatta che si sarebbe ottenuta in assenza di errori e x + δx la soluzionefornita dal calcolatore. Generalmente la matrice A non si conosce esattamente perche icoefficienti possono essere delle misure, oppure sono calcolati con errori di arrotondamento,quindi nella pratica ci troviamo a lavorare con un sistema del tipo

Ax = b

ove δA = A − A. E naturale chiedersi quanto x e vicino ad x. Questo e un problema diteoria delle perturbazioni matriciale.

Innanzi tutto occorre trovare le condizioni sotto le quali la matrice A risulta non singolare.

TeoremaSia A non singolare. Se

‖A−1δA‖ < 1

oppure

‖A−1‖‖δA‖ < 1

allora A + δA e non singolare

36


Per semplita consideriamo inizialmente il caso in cui sia perturbato il solo termine noto.Sia δb ∈ Rn il vettore che perturba il termine noto. Allora la soluzione sara x + δx, piuprecisamente:

A(x + δx) = b + δb

da cui si ricavaAδx = δb, δx = A

−1δb

Per una qualsiasi norma matriciale indotta dalla norma vettoriale si ottiene:

‖δx‖ = ‖A−1δb‖ ≤ ‖A−1‖‖δb‖

d’altra parte abbiamo‖b‖ = ‖Ax‖ ≤ ‖A‖‖x‖

combinando insieme queste ultime due relazioni si ottiene:

‖δx‖‖x‖ ≤ ‖A‖‖A−1‖‖δb‖

‖b‖

Generalmente si indica con µ(A) = ‖A‖‖A−1‖ il numero di condizionamento della matriceA, si osservi che

µ(A) = ‖A‖‖A−1‖ ≥ µ(A) = ‖AA−1‖ = 1

quindi se µ(A) assume valori piccoli allora a piccole perturbazioni sui dati corrispondonopiccole perturbazioni sulla soluzione e quindi il problema e ben condizionato, cioe la matricedel sistema e ben condizionata. Se µ(A) assume valori grandi allora piccole variazioni sui datipossono indurre grandi perturbazioni nella soluzione, e quindi il problema e mal condizionato.

Nel caso in cui consideriamo le perturbazioni anche sui coefficienti della matrice Asupposto che la matrice A + δA sia non singolare, indicata ancora con x + δx la soluzionedel sistema perturbato

(A + δA)(x + δx) = b + δb

37


risulta:

‖δx‖‖x‖ ≤ µ(A)

‖δA‖‖A‖ +

‖δb‖‖b‖

1− µ(A)‖δA‖‖A‖

se indichiamo con εx la perturbazione relativa sulla soluzione indotta dalle perturbazionirelative dei datie εA, εb allora

εx ≤ µ(A)εA + εb

1− µ(A)εA

Se consideriamo la matrice di Hilbert H i cui elementi sono dati da:

hi,j =1

i + j − 1, i, j = 1, . . . , n

La matrice inversa esatta ( ottenuta per via analitica) e data da

h−1i,j

(−1)i+j(n + i− 1)!(n + j − 1)!

(i + j − 1)[(i− 1)!(j − 1)!]2(n− i)!(n− j)!

Il numero di condizione e stimato

µ2(H) = O(eαn

)

con α costante maggiore di zero.

38


x

x^

x

x^

39


Metodi Iterativi per Sistemi Lineari

Per risolvere un sistema n × n lineare Ax = b si possono utilizzare anche i metodiiterativi, che risultano particolarmente convenienti se la matrice A e sparsa, cioe se il numerodegli elementti non nulli di A e dell’ordine di n. Se si usa un metodo diretto per matricisparse si ha il cosı detto fill in cioe elementi della matrice, che inizialmente sono nulli, nelcorso della fattorizzazione diventano 6= 0. Quindi se A e sparsa, di grandi dimensioni e nonha struttura allora i metodi iterativi risultano convenienti.

DefinizioneUna successione {x(k)} di vettori in Rn si dice convergente al vettore x(∗) di Rn seesiste una norma per cui risulta

limk→∞

‖x(k) − x(∗)‖ = 0

in tal caso si pone limk→∞ x(k) = x(∗). La condizione di convergenza si traduce inuna condizione di convergenza delle successioni formate dalle singole componenti. Infatticonsiderando la norma ∞, poiche

|x(k)i

− x(∗)i| ≤ ‖x(k) − x

(∗)‖, i = 1, . . . , n,

da cui si ha limk→∞ |x(k)i

− x(∗)i| = 0 e quindi

limk→∞

x(k)i

= x(∗)i

, i = 1, . . . , n,

vale inoltre il viceversa.

40


Per le successioni di matrici {A(k)} si puo dare una definizione di convergenza analoga.

Di fondamentale importanza per la convergenza dei metodi iterativi sono i seguentiteorema:

TeoremaSia A ∈ Rn×n allora

limk→∞

Ak

= 0 se e solo se ρ(A) < 1

dove ρ(A) e il raggio spettrale della matrice A

TeoremaSia A ∈ Rn×n allora

det(I − A) 6= 0 e limk→∞

kX

i=0

Ai

= (I − A)−1

se e solo seρ(A) < 1

41


Generalita sui metodi iterativi

Sia A ∈ Rn×n non singolare e si consideri la decomposizione di A nella forma

A = M −N

dove M e una matrice non singolare. Sostituendo nel sistema lineare

Mx−Nx = b

risultax = M

−1Nx + M

−1b

PostoP = M

−1N e q = M

−1b,

si ottiene il seguente sistemax = Px + q

equivalente al sistema originale.

Dato un vettore iniziale x(0) si considera la successione x(1), x(2), . . . , cosı definita

x(k)

= Px(k−1)

+ q, k = 1, 2, . . .

Se la successione x(k) e convergente e si indica con

x(∗)

= limk→∞

x(k)

,

42


allora passando al limite risulta

x(∗)

= Px(∗)

+ q

cioe x(∗) e la soluzione del sistema.La matrice P si dice matrice di iterazione del metodo.

Al variare del vettore iniziale x(0) si ottengono diverse successioni {x(k)}, alcune dellequali possono essere convergenti ed altre no.

Un metodo iterativo e detto convergente se, qualunque sia il vettore iniziale x(0), la

successione {x(k)} e convergente.

TeoremaIl metodo iterativo e convergente se e solo se ρ(P ) < 1.Infatti

x(∗) − x

(k)= P (x

(∗) − x(k−1)

), k = 1, 2, . . .

indicato con

e(k)

= x(∗) − x

(k),

il vettore errore alla k-esima iterazione, si ha

e(k)

= Pe(k−1)

, k = 1, 2, . . .

e quindi

e(k)

= Pe(k−1)

= P2(e

(k−2)) = . . . = P

ke(0)

.

Se ρ(P ) < 1 allora

limk→∞

P(k)

= 0,

e quindi

limk→∞

e(k)

= 0.

Viceversa e facile verificare che se l’autovalore massimo di A e minore di 1 allora l’erroree(k) → 0 per k tendente all’infinito.

43


Se esiste una norma matriciale indotta ‖ · ‖ per cui ‖P‖ < 1 il metodo iterativo econvergente.

Infatti

ρ(P ) ≤ ‖P‖sia λ un autovalore e x il corrispondente autovettore normalizzato rispetto alla norma ‖ · ‖,

Ax = λx, ‖x‖ = 1

allora passando alle norme

|λ| = ‖Ax‖,da cui segue

|λ| ≤ max‖v‖=1

‖Av‖ = ‖A‖

questa relazione vale per ogni autovalore λ di A e quindi anche per quello di modulo massimo.

Osservazioni

• Se il determinate della matrice di iterazione P e ≥ 1 allora almeno uno degli autovaloridi P e di modulo ≥ 1 e quindi il metodo iterativo non e convergente.

• poiche la traccia e uguale alla somma degli autovalori, allora se |tr(P )| ≥ n almeno unautovalore di P e di modulo maggiore o uguale a 1 e quindi il metodo non e convergente(si ricorda anche che tr(P ) =

Pni=1 pi,i).

Quindi condizioni necessarie affinche il metodo iterativo sia convergente:

• |det(P )| < 1

• |tr(P )| < n

44


Controllo della Convergenza

Analizziamo l’errore da cui e affetto il vettore x(k) rispetto alla soluzione x(∗) soluzione delsistema:

‖e(k)‖ ≤ ‖Pk‖‖e(0)‖,quindi ‖Pk‖ esprime la riduzione, rispetto all’errore iniziale, dell’errore al k-esimo passo.Questa misura risulta inadatta per una valutazione della velocita di convergenza di un metodoche sia indipendente dal numero delle iterazioni. Infatti se P e Q sono due matrici di iterazioneassociate a 2 metodi diversi, puo accadere che per una norma ‖ · ‖ particolare esistano dueindici interi j 6= k tali che

‖Pk‖ < ‖Qk‖ ‖Pj‖ > ‖Qj‖

Esempio

P =

»0.5 00 0.6

–, Q =

»0.5 0.250 0.5

–

si ha

Pk

=

"0.5k 0

0 0.6k

#, Q

k=

"0.5k k0.25k

0 0.5k

#

Utilizzando la norma ∞ risulta

‖Pk‖∞ = 0.6k

e ‖Qk‖∞ = (2 + k)0.5k

per k = 1 : 9 si ha ‖Pk‖ < ‖Qk‖ e per k ≥ 10 ‖Pk‖ > ‖Qk‖, se si usa la normadue il cambiamento si ha per k = 7.

Quindi si percorre un’altra strada. Se e(k−1) 6= 0, la quantita

‖e(k)‖‖e(k−1)‖

45


esprime la riduzione dell’errore al k-esimo passo.La media geometrica delle riduzioni dell’errore sui primi k − passi

σk = k

vuut‖e(1)‖‖e(0)‖

‖e(2)‖‖e(1)‖

. . .‖e(k)‖‖e(k−1)‖

≤

k

vuut‖e(k)‖‖e(0)‖

esprime la riduzione media per passo dell’errore relativo ai primi k passi, quindi dalla relazioneprecedente ricaviamo che

σk ≤kq‖Pk‖

dove il segno di uguale vale solo per particolari valori di e(0). Se facciamo tendere k all’infinitosi ottiene la riduzione asintotica media per passo che come dimostra il seguente teorema nondipende dalla norma usata.

TeoremaSia A ∈ Rn×n e sia ‖ · ‖ una norma indotta. Allora

limk→∞

kq‖P (k)‖ = ρ(P )

La quantita ρ(P ) e indipendente dalla norma e dall’indice di iterazione k, viene quindi assuntacome misura della velocita di convergenza del metodo.Il numero k di iterazioni richieste per ridurre l’errore di 1

10 (cioe approssimativamente perottenere una cifra decimale in piu) e

[ρ(P )]k ∼ 1

10, da cui k ∼ − 1

log10(ρ(P ))

46


Si definisce tasso asintotico di convergenza del metodo iterativo la costante

R = −log10(ρ(P ))

Puo essere utile la localizzazione degli autovalori di una matrice, che si puo otteneremediante l’analisi dei cerchi di Gershgorin.

47


Criteri di Arresto

In generale con un metodo iterativo non e possibile calcolare la soluzione in un numerofinito di iterazioni, occorre quindi determinare dei criteri per l’arresto del procedimento. Sonoanaloghi a quelli visti nella ricerca degli zeri di funzioni non lineari.Fissata una tolleranza ε che tiene conto anche della precisione utilizzata nei calcoli siconsiderano i seguenti test di arresto:

‖x(k) − x(k−1)‖ ≤ εa,

se x(k) 6= 0

‖x(k) − x(k−1)‖‖x(k)‖

≤ εr

Queste condizioni non garantiscono che la soluzione sia stata approssimata con precisione ε.Per vedere cio ricordiamo il seguente risultato:se per una qualche norma matriciale indotta risulta ‖A‖ < 1

⇒ ∃ (I + A)−1

e

‖(I + A)−1‖ ≤ 1

1− ‖A‖Allora consideriamo

x(k+1) − x

(k)= [x

(∗) − x(k)

]− [x(∗) − x

(k+1)] =

e(k) − e

(k+1)= (I − P )e

(k)= (I + (−P ))e

(k)

48


cioe

e(k)

= (I + (−P ))−1

(x(k+1) − x

(k))

e passando alle norme

‖e(k)‖ ≤ ‖(I + (−P ))−1‖‖x(k+1) − x

(k)‖ ≤

≤ 1

1− ‖P‖‖x(k+1) − x

(k)‖

questo significa che il test di arresto

‖x(k+1) − x(k)‖ ≤ εa + εr‖x(k)‖

puo essere verificato, ma l’errore ‖e(k)‖ puo essre molto grande. Questo accade quando‖P‖ ∼ 1 per cui

1

1− ‖P‖ ∼ ∞.

In un programma che implementa un metodo iterativo deve essere comunque previsto uncontrollo per interrompere l’esecuzione quando il numero delle iterazioni diventa troppoelevato.Puo accadere che un metodo iterativo la cui matrice di iterazione P e tale che ρ(P ) < 1,per gli effetti degli errori di arrotondamento nella pratica non converga. Questo accade inparticolare quando la matrice A e fortemente mal condizionata e ρ(P ) e vicino a 1.Un primo stimatore e costituito dal residuo che ad ogni iterazione e definito da

rk

= b− Axk

Ponendo x = xkmin per il quale si abbia

‖b− Axkmin‖ ≤ ε‖b‖

49


si ha

‖e(kmin)‖‖x(∗)‖

≤ εµ(A)

Il controllo del residuo e significativo solo se il numero di condizione di A e ragionevole.Un metodo iterativo rispetto ad un metodo diretto e in generale meno sensibile alla

propagazione degli errori. Infatti il vettore x(k) puo essere considerato come il vettore

generato da una sola iterazione a partire dal vettore iniziale x(k−1) e quindi risulta affettodagli errori di arrotondamento generati dalla sola ultima iterazione.

In un metodo iterativo ad ogni iterazione il costo computazionale e principalmentedeterminato dalla operazione di moltiplicazione della matrice P per un vettore, che richiede

n2 operazioni moltiplicative se la matrice A non ha specifiche proprieta.Se A e sparsa, cioe ha un numero di elementi non nulli dell’ordine di n, la moltiplicazione diP per un vettore richiede un numero di operazioni moltiplicative dell’ordine di n.In questo caso i metodi iterativi possono risultare competitivi con quelli diretti.

Metodo Iterativo di Jacobi

E anche noto come il metodo delle sostituzione simultanee, oppure di Richardson.

Corrisponde a una particolare decomposizione della matrice dei coefficienti A ∈ Rn×n,

det(A) 6= 0, nella somma di due matrici M , N ∈ Rn×n, det(M) 6= 0,

A = M −N, M = diag(A), N = {−ai,j}i6=j

dunque

M =

26664

a1,1 0 . . . 0

0 a2,2 . . . 0. . .

0 . . . 0 an,n

37775

N =

26664

0 −a1,2 . . . −a1,n−a2,1 0 . . . −a2,n

. . .

−an,1 . . . −an,n−1 0

37775

50


La matrice di iterazione PJ del metodo di Jacobi e allora

PJ =

266666664

0 −a1,2a1,1

. . . −a1,na1,1

−a2,1a2,2

0 . . . −a2,na2,2

. . .

−an,1an,n . . . −an,n−1

an,n 0

377777775

da cui il metodo di Jacobi

x(k+1)

= PJx(k)

+ q, q = [b1

a1,1,

b2a2,2

, . . . ,bn

an,n]T

ossia

x(k+1)i

=1

ai,i

264bi −

nX

j=1,j 6=i

ai,jxj(k)

375

l’implementazione del metodo richiede due vettori, xold, xnew; alla fine di ogni ciclo i se i

test non sono verificati, si aggiorna xnew → xold, in modo da poter riutilizzare i due vettoriper il passo successivo. Le singole componenti xnew

i del vettore xnew vengono costruite inmodo indipendente tra loro, il che rende parallelo l’algoritmo.

51


%JACOBI% x=Jacobi(A,b,x0,epsilon)% Pre: A non singolare, A(i,i) non nulli per i=1:n.% In realta’ perche’il problema sia ben condizionato,% A deve% essere a diagonale dominante.% La funzione restituisce in x la soluzione del sistema Ax=b,% risolto con il metodo semiiterativo di Jacobi.% Riceve in x0 il vettore di innesco e in epsilon la tolleranza% desiderata.%% See also GAUSS_SEIDELfunction x=Jacobi(A,b,x0,maxiter,psilon) n=length(b); x=x0;iter=0;r=b-A*x;

while norm(r)>epsilon*norm(b)&iter<maxiteriter=iter+1;

for i=1:nsom=0;for j=1:n

if i~=jsom=som+A(i,j)*x0(j);

endendsom=(b(i)-som)/A(i,i);x(i)=som;

endx0=x;r=b-A*x;

end return

L’algoritmo di Jacobi e definito solo se ai,i 6= 0 ∀i = 1 : n, e converge se e solo se

ρ(PJ ) < 1.Se tale ipotesi non e verificata, ma la matrice A e non singolare, e possibile riordinare leequazioni e le incognite del sistema Ax = b in modo da ottenere un sistema equivalente alsistema originale

Ax = b

52


per il quale il metodo di Jacobi risulti definito. Occorre pero che le condizioni di convergenzasiano verficate.

Esempio

Applichiamo il metodo di Jacobi al sistema lineare Ax = b

A =

24

20 2 −12 13 −21 1 1

35 , x =

24

x1x2x3

35 , b =

24

25302

35 .

La decomposizione di Jacobi e

M =

24

20 0 00 13 00 0 1

35 , N =

24

0 −2 1−2 0 2−1 −1 0

35 .

La matrice ed il vettore di iterazione di Jacobi sono

PJ =

264

0 − 110

120

− 213 0 2

13−1 −1 0

375 , q =

264

5430132

375 .

Gli autovalori di PJ sono λ1 = 0.114487, λ2,3 = −0.0572433 ± 0.4453i allora

ρ(PJ ) = |λ2,3| = 0.44896 < 1. Il metodo di Jacobi converge alla soluzione

x = [1, 2.− 1]T .

Permutiamo ora le prime due equazioni e consideriamo il sistema equivalente Ax = b:

A =

24

2 13 −220 2 −11 1 1

35 , x =

24

x1x2x3

35 , b =

24

30252

35 .

53


Otteniamo

M =

24

2 0 00 2 00 0 1

35 , N =

24

0 −13 2−20 0 1−1 −1 0

35 .

PJ =

264

0 −132 1

−10 0 12

−1 −1 0

375 , q =

24

15252

2

35 .

Gli autovalori della matrice PJ sono:

λ1 = 8.07104 + 8.88178× 10−16i,

λ2 = −0.208805− 1.77636× 10−15i e

λ3 = −7.86223− 8.88178× 10−16ied in valore assoluto|λ1| = 8.07104 allora ρ(PJ ) = |λ1| > 1. Il metodo di Jacobi in questo caso nonconverge.

54


Metodo iterativo di Gauss-Seidel

E anche noto come metodo delle sostizioni successive, oppure di Liebmann.

Corrisponde alla seguente decomposizione della matrice dei coefficienti A ∈ Rn×n,

det(A) 6= 0, nella somma di due matrici M , N ∈ Rn×n, det(M) 6= 0,

A = M −N, M = diag(A) + {ai,j}i>j

N = {−ai,j}i<j = −UA

quindi

M =

26664

a1,1 0 . . . 0

0 a2,2 . . . 0. . .

0 . . . 0 an,n

37775+

26664

0 0 . . . 0a2,1 0 . . . 0

. . . . . .

an,1 . . . an,n−1 0

37775 =

=

26664

a1,1 0 . . . 0

a2,1 a2,2 . . . 0. . . . . .

an,1 . . . an,n−1 an,n

37775

e

N =

266664

0 −a1,2 . . . −a1,n0 0 . . . −a2,n

. . . . . .

0. . . 0 0

377775

La matrice ed il vettore di iterazione del metodo non sono immediatamente esprimibili in formaesplicita sono dati da

PGS = M−1

N, q = M−1

b

55


E possibile invece esplicitare

x(k+1)

= PGSx(k)

+ q

Mx(k+1)

= Nx(k)

+ b

DAx(k+1)

+ LAx(k+1)

= −UAx(k)

+ b

DAx(k+1)

= −LAx(k+1) − UAx

(k)+ b

x(k+1)

= −D−1A LAx

(k+1) −D−1A UAx

(k)+ D

−1A b

D−1A LA =

2666664

0 0 . . . 0a2,1a2,2

0 . . . 0

. . .an,1an,n . . .

an,n−1an,n 0

3777775

,

D−1A UA =

26666664

0a1,2a1,1

. . .a1,na1,1

0 0 . . .a2,na2,2

. . .

0 0 . . . 0

37777775

ed analogamente per

D−1A b = [

b1a1,1

,b2

a2,2, . . . ,

bn

an,n]T

ossia

x(k+1)i

=1

ai,i{bi −

i−1X

j=1

ai,jx(k+1)j

−i+1X

j=1

ai,jx(k)j}

, i = 1, . . . , n

L’implementazione del metodo di Gauss-Seidel richiede un solo vettore xcurr, le singolecomponenti xcurr

i vengono utilizzate non appena sono state calcolate, per cui in questo casol’algoritmo e sequenziale.

56


In Matlab lo si puo implementare come segue:

%GAUSS_SEIDEL% function [x,count]=Gauss_Seidel(A,b,x0,epsilon,upper)% Pre: A non singolare. Perche’ il problema sia ben% condizionato dev’essere A a diagonale dominante%% La function risolve il sistema lineare Ax=b con il metodo% di Gauss Seidel.% Riceve in x0 il vettore di innesco, in epsilon l’accuratezza% desiderata, in upper un limite superiore sul numero di% passi da eseguire.% Restituisce in x la soluzione del sistema e in count il numero% di passi che sono stati effettivamente necessari per arrivare% alla soluzione. Se passi=upper si consiglia di riprovare con% un valore di upper piu’ alto.%%See also DD, JACOBI

function [x,count]=GaussSeidel(A,b,x0,epsilon,upper) n=length(b);x=x0; r=b-A*x; count=0; while norm(r)>=epsilon*norm(b);

x_prec=x;for i=1:n

x(i)=(b(i)-A(i,i+1:n)*x(i+1:n)-A(i,1:i-1)*x(1:i-1))/A(i,i);endcount=count+1;if ((norm(x-x_prec)<epsilon) & (count>=upper))

return;endr=b-A*x;

end return;

L’algoritmo di Gauss-Seidel e definito solo se ai,i 6= 0 ∀i = 1 : n, e converge se e

solo se ρ(PGS) < 1.Se tale ipotesi non e verificata, ma la matrice A e non singolare, e sempre possibile riordinarele equazioni del sistema lineare in modo da ottenere un sistema equivalente in cui Gauss-Seidel

57


sia definito.

EsempioApplichiamo il metodo di Gauss-Seidel al sistema Ax = b

A =

24

20 2 −12 13 −21 1 1

35 , x =

24

x1x2x3

35 , b =

24

25302

35 .

La decomposizione di Gauss-Seidel e

M =

24

20 0 02 13 01 1 1

35 , N =

24

0 −2 10 0 20 0 0

35 .

La matrice ed il vettore di iterazione di Gauss-Seidel sono

PGS =

264

0 − 110

120

0 165

9130

0 11130 − 51

260

375 , q =

264

54

5526

−7152

375 .

Gli autovalori di PGS sono λ1 = 0.0λ2 = −0.243858λ3 = 0.0630885allora ρ(PGS) = |λ2| < 1, il metodo di Gauss-Seidel converge alla soluzione x =

[1, 2,−1]T .In questo esempio vale

ρ(PGS) = 0.243858 < ρ(PJ ) = 0.44896 < 1

quindi per questa matrice il procedimento di Gauss-Seidel converge piu velocemente di quellodi Jacobi. Questo spesso e vero, ma non sempre.

58


• Esistono matrici per le quali il metodo di Gauss-Seidel converge, mentre Jacobi diverge eviceversa.

• Esistono matrici per le quali il metodo di Gauss-Seidel converge piu velocemente di quellodi Jacobi e viceversa.

• non e lecito utilizzare la norma della matrice di iterazione per decidere quale fra idue metodi converge piu velocemente, si deve utilizzare il raggio spettrale, perche puoaccadere che

ρ(PGS) < ρ(PJ ) & ‖PJ‖∞ ≤ ‖PGS‖∞

oppure viceversa.

Convergenza dei metodi di Jacobi e di Gauss-Seidel

Le condizioni necessarie e sufficienti ρ(P ) < 1 oppure solo sufficienti ‖P‖ < 1 vistein precedenza valgono anche per questi due metodi. In particolare per questi due metodiesistono altre condizioni (necessarie e sufficienti o solo sufficienti) per la convergenza basatisulle proprieta della matrice dei coefficienti A, in genere facilmente verificabili.

• Matrici a Diagonale Dominante Stretta;

• Matrici Irriducibili a Diagonale Dominante

• Matrici Hermitiane definite positive

Per matrici con diagonale reale non nulla e con Pj non negativa (cioe pi,j ≥ 0 ∀i, j) alloravale il seguente teorema:

59


Teorema di Stein RosenbergSia A una matrice con elementi diagonali non nulli; inoltre sia Pj non negativa. Allora valeuna sola delle relazioni seguenti:

• 1. ρ(PGS) = ρ(PJ ) = 0

• 2. ρ(PGS) < ρ(PJ ) < 1

• 3. ρ(PGS) = ρ(PJ ) = 1

• 4. ρ(PGS) > ρ(PJ ) > 1

Pertanto, per una matrice A che soddisfa le ipotesi del teorema per studiare la convergenza(o divergenza) del metodo di Gauss-Seidel e sufficiente studiare la convergenza (o divergenza)del metodo di Jacobi (piu semplice). Nel caso di convergenza sappiamo che il metodo diGauss-Seidel converge piu velocemente di quello di Jacobi.

Matrici TridiagonaliUna matrice si dice tridiagonale se ai,j = 0, per tutti gli indici i, j tali che |i− j| > 1.

Teorema

Sia A una matrice tridiagonale, con elementi diagonali non nulli. Allora valgono leseguenti relazioni

• 1. se λ e autovalore di PJ allora λ2 e autovalore di PGS ;

• 2. se λ e autovalore di PGS allora√

λ e autovalore di PJ ;

• 3. ρ(PGS) = ρ(PJ )2. Quindi per una matrice tridiagonale con diagonale non nulla, ilmetodo di Gauss-Seidel ha una velocita di convergenza asintotica doppia rispetto a quelladel metodo di Jacobi.

Riassumendo

• 1. Un metodo iterativo per la soluzione di un sistema lineare costruisce, a partire

da un vettore iniziale x(0), una successione di vettori x(k) sui quali si richiede che

limk→∞ x(k) = x(∗);

60


• 2. Condizione necessaria e sufficiente affinche un metodo iterativo converga per ogni

possibile scelta di x(0) e che il raggio spettrale della matrice di iterazione sia minore di1.

• 3. I metodi iterativi piu classici sono quelli di Jacobi e Gauss-Seidel. Condizione sufficienteper la loro convergenza e che la matrice di iterazione sia a dominanza diagonale stretta osimmetrica definita positiva.

• 4. tre sono i possibili criteri di arresto, si possono compattare i due soli, il primo basatosul controllo del residuo e significativo se la matrice e ben condizionata, il secondo basatosul numero max di ietrazioni.Essi vanno tutti verificati ad ogni iterazione.

Hermitiane

definitepositive

stretta

Gauss−Seidel

convergente

Predominanza diagonale

Per quanto riguarda i metodi iterativi esiste un’ampia famiglia di metodi moderni, noticome metodi di Krylov, decisamente piu efficienti di quelli qui proposti. Tra essi spiccano perle proprieta di terminazione finita, cioe per raggiungere in aritmetica esatta la soluzione esattadel sistema in un numero finito di iterazioni, i metodi del gradiente coniugato e del GMRES

61


(Generalized Minimum RESidual); il primo e applicabile ai soli sistemi con matrice simmetricadefinita positiva. Entrambi i metodi sono implememntati in Matlab all’interno del toolboxsparfun e sono richiamabili con i comandipcggmresUn metodo iterativo, per quanto efficiente, convergera lentamente se la matrice dei coefficientie mal condizionata. Sono state implementate numerose startegie di precondizionamento,alcune delle quali puramente algebriche basate sulla fattorizzazione LU incomplete. Sonoimplementate in Matlab nella funzioneluinc

ESERCIZIOSi consideri il sistema lineare la cui matrice dei coefficienti e data da:

A =

» −10 2β 5

–

Si dia una condizione sufficiente sul numero reale β affinche i metodi di Jacobi e di Gauss -

Seidel convergano entrambi.

62

Risoluzione di Sistemi lineari - dm.unibo.itguerrini/html/an_08_09/sistemi_lineari_completi.pdf ·...

Documents

Transcript of Risoluzione di Sistemi lineari - dm.unibo.itguerrini/html/an_08_09/sistemi_lineari_completi.pdf ·...