Il metodo di Eliminazione di Gauss con la tecnica del pivoting

Pivoting e stabilitaLa fattorizzazione LU

Metodo di CholeskyLe function Matlab

Corso di Matematica per la Chimica

Dott.ssa Maria Carmela De Bonisa.a. 2013-14

Dott.ssa Maria Carmela De Bonis a.a. 2013-14 Corso di Matematica per la Chimica



Pivoting e stabilita

Se la matrice A non appartiene a nessuna delle categorie precedentipuo accadere che al k−esimo passo risulti a(k)k,k = 0, e quindi ilmetodo di Gauss, cosı com’e stato descritto, non puo proseguire.

Tuttavia se a(k)k,k = 0 necessariamente qualche altro elemento

a(k)i,k , i = k + 1, . . . , n, della colonna k−esima della matrice deicoefficienti deve essere non nullo, altrimenti la matrice dei coefficientisarebbe singolare.

Se, ad esempio, a(k)r,k 6= 0, basta scambiare l’equazione k−esima con lar−esima e poi procedere con le eliminazioni.

Dunque ogni sistema non singolare, mediante opportuni scambi dirighe, puo essere sempre ricondotto alla forma triangolare superiorecon il metodo di Gauss.




Discutiamo ora della stabilita del metodo di Gauss.

E possibile provare che l’esistenza di pivot molto piccoli (in valoreassoluto) rispetto all’ordine di grandezza degli elementi della matricee causa di cancellazione numerica e quindi di instabilita.

Dunque per assicurare una migliore stabilita numerica al metodo dieliminazione di Gauss e possibile permutare l’ordine delle equazionianche quando l’elemento pivot non e nullo ma e “piccolo”.

Tale strategia algoritmica e detta pivoting.




Il pivoting consiste nello scegliere, al generico passo k−esimo,l’elemento pivot in maniera ottimale.

Le due strategie di pivoting piu utilizzate sono le seguenti:

pivoting parziale: si sceglie r uguale al piu piccolo intero ≥ ktale che

|a(k)r,k | = maxk≤i≤n

|a(k)i,k |

e, se r 6= k, si scambia l’equazione k−esima con l’r−esima;

pivoting totale: si sceglie la coppia (r, s), con r, s ≥ k tale che

|a(k)r,s | = maxk≤i,j≤n

|a(k)i,j |

e si scambiamo l’equazione k−esima con l’r−esima el’incognita k−esima (con il suo coefficiente) con l’s−esima.




La strategia del pivoting totale combinata con il metodo di Gaussassicura la stabilita dell’algoritmo complessivo.

Tuttavia essa puo risultare molto costosa.

La strategia di pivoting parziale e meno “costosa” e, poiche, ingenerale, risulta soddisfacente nella maggior parte dei casi, essa e lastrategia piu utilizzata.

Osservazione

Il metodo di eliminazione di Gauss senza pivoting e comunquenumericamente stabile quando:

la matrice A del sistema e a diagonale dominante per colonne (inquesto caso la strategia di pivot non produce scambi);

la matrice A del sistema e simmetrica e definita positiva (inquesto caso la strategia effettua scambi ma non producemiglioramenti).




Esempio

Consideriamo il sistema lineare Ax = b di ordine n = 18, dove

ai,j = cos((j− 1)

2i− 12n

π

), i, j = 1, . . . , n,

e

bi =

n∑j=1

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

la cui soluzione esatta e x = (1, 1, . . . , 1)T .

Tale matrice e ben condizionata, risultando

cond(A) = ‖A‖∞‖A−1‖∞ = 16.90251471518910

Risolvendo il sistema con il metodo di eliminazione di Gauss senza econ la strategia di pivoting parziale, si ottengono i seguenti risultati:




Gauss Gauss + pivoting parziale9.999999968425205e-001 9.999999999999993e-0011.000000005848771e+000 1.000000000000000e+0009.999999953310822e-001 1.000000000000001e+0001.000000003291504e+000 1.000000000000001e+0009.999999977773912e-001 9.999999999999997e-0011.000000001703849e+000 1.000000000000000e+0009.999999983501708e-001 1.000000000000000e+0001.000000001767299e+000 1.000000000000001e+0009.999999982508074e-001 9.999999999999998e-0011.000000001433063e+000 9.999999999999999e-0019.999999991425489e-001 9.999999999999993e-0011.000000000200525e+000 1.000000000000000e+0001.000000000347843e+000 9.999999999999996e-0019.999999993163441e-001 9.999999999999990e-0011.000000000800448e+000 9.999999999999996e-0019.999999992564082e-001 1.000000000000000e+0001.000000000564552e+000 1.000000000000000e+0009.999999996959291e-001 9.999999999999999e-001




La fattorizzazione LU

Il metodo di eliminazione di Gauss, dal punto di vista matriciale, puoessere riletto come la costruzione di una successione di matrici

[A|b] = [A(1)|b(1)], . . . , [A(k)|b(k)], . . . , [A(n)|b(n)]

in modo tale che A(n) sia triangolare superiore e b(n) sia il nuovotermine noto.

Le matrici della successione sono tra loro legate da unatrasformazione del tipo

[A(k+1)|b(k+1)] = M(k)[A(k)|b(k)], k = 1, . . . , n− 1




dove

M(k) =

1 0 . . . . . . 0

0 1 0...

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

0 0 1...

0 0 mk+1,k 1. . .

0 0 mk+2,k 0. . .

......

......

. . . . . . 00 . . . 0 mn,k 0 . . . 0 1

e detta matrice elementare di Gauss.




Si ha quindi che

A = [M(1)]−1A(2) = [M(1)]−1[M(2)]−1A(3) = · · ·= [M(1)]−1 · · · [M(n−1)]−1A(n)

eb = [M(1)]−1 · · · [M(n−1)]−1b(n).

Ponendo

L = [M(1)]−1 · · · [M(n−1)]−1, U = A(n) e y = b(n),

si ottieneA = LU e b = Ly,

dove U e la matrice triangolare superiore che si ottiene alla fine delmetodo di eliminazione di Gauss e L e una matrice triangolareinferiore (e prodotto di inverse di matrici triangolari inferiori).




La matrice L e definita come

L =

1 0 . . . . . . 0

−m2,1 1 0...

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

. . . 1...

... −mk+1,k 1. . .

−mk+2,k. . . . . .

......

. . . 1 0−mn,1 . . . . . . −mn,k . . . . . . −mn,n−1 1

e dunque per costruirla basta, ad ogni passo del metodo di Gauss,memorizzare i moltiplicatori cambiati di segno.




schema algoritmo

for k=1:n-1for i=k+1:n

ai,k = ai,k/ak,k; costruzione matrice Lfor j=k+1:n

ai,j = ai,j − ai,k ∗ ak,j; costruzione matrice Uend

endendL=eye(n)+tril(A,-1)U=triu(A)




Dunque, con piccoli accorgimenti algoritmici, dal metodo di Gauss epossibile calcolare le due matrici L ed U tali che

A = LU.

Si effettua in tal modo una decomposizione della matrice A dettafattorizzazione LU di A.




Calcolate le matrici L e U, il sistema Ax = b, puo essere risoltomediante i due sistemi

Ly = b e Ux = y,

il primo triangolare inferiore e il secondo triangolare superiore.

Risolvere il primo sistema equivale a calcolare il nuovo termine noto.La differenza con l’algoritmo di Gauss e che, invece di farlocontestualmente alla riduzione in forma triangolare, tale calcolo vieneeffettuato in modo indipendente.

Il costo computazionale complessivo ammonta dunque a n3

3 + n2 dicui n3

3 per la fattorizzazione LU e n2 per le due sostituzioni(all’indietro e in avanti).




Le fattorizzazioni di matrici hanno diversi utilizzi.

Nel caso specifico della fattorizzazione LU la prima immediataapplicazione e il calcolo del determinante di A.

Infatti

det(A) = det(LU) = det(L) det(U) = det(U) =

n∏i=1

ui,i.




Un’altra possibile applicazione scaturisce dall’esigenza di risolvere psistemi che hanno tutti la stessa matrice dei coefficienti A, cioe

Ax1 = b1, Ax2 = b2, . . . , Axp = bp,

o in altri termini il sistema

AX = B, con A ∈ Rn×n, X,B ∈ Rn×p.

Infatti in tal caso la fattorizzazione viene effettuata una sola volta eogni sistema viene poi risolto mediante 2 algoritmi di sostituzione,per un costo computazionale complessivo che si riduce quindi an3

3 + pn2 (contro p(

n3

3 + n2

2

)se si risolvesse ciascun sistema

indipendentemente dagli altri).

Osserviamo infine che se p = n e B = I, risolvere il sistema AX = B,e equivalente a calcolare A−1.




Il metodo di Gauss con la variante del pivoting esegue ancora unafattorizzazione di matrice nei seguenti termini

PA = LU e Pb = Ly

dove P ∈ Rn×n, detta matrice di permutazione, contiene leinformazioni relative agli scambi di righe.

Vale il seguente

TeoremaPer ogni matrice A ∈ Rn×n esiste una matrice di permutazioneP ∈ Rn×n tale che

PA = LU.




Metodo di Cholesky

Sia A ∈ Rn×n una matrice simmetrica definita positiva.

In questo caso speciale e possibile costruire un algoritmo, antagonistadel metodo di Gauss, e quindi una fattorizzazione alternativa a quellaLU, che e piu conveniente dal punto di vista computazionale.

Vale il seguente

TeoremaSe A ∈ Rn×n e una matrice simmetrica definita positiva esiste ed eunica la fattorizzazione

A = LLT

dove L e una matrice triangolare inferiore con elementi diagonali nonnulli.




Posto A = (ai,j)i,j=1,...,n e L = (li,j)i,j=1,...,n, si ha

ai,j =

n∑k=1

li,klTk,j =n∑

k=1

li,klj,k.

Poiche la matrice e simmetrica, possiamo considerare solo glielementi di A con j ≤ i. Otteniamo

ai,j =

j−1∑k=1

li,klj,k + li,jlj,j, i = 1, . . . , n, j = 1, . . . , i− 1

e

ai,i =

i−1∑k=1

l2i,k + l2i,i, i = 1, . . . , n.




Da cuil1,1 =

√a1,1

li,j =1lj,j

[ai,j −

j−1∑k=1

li,klj,k

]i = 2, . . . , nj = 1, . . . , i− 1

li,i =

√√√√ai,i −i−1∑k=1

l2i,k




schema algoritmo

l1,1 =√a1,1

for i=2:nfor j=1:i-1

li,j = 0;for k=1:j-1

li,j = li,j + li,k ∗ lj,k; j− 1 operazioniendli,j = (ai,j − li,j)/lj,j; 1 operazione

endli,i = 0;for k=1:i-1

li,i = li,i + l2i,k; i− 1 operazioniendli,i =

√ai,i − li,i; 1 operazione

endDott.ssa Maria Carmela De Bonis a.a. 2013-14 Corso di Matematica per la Chimica



costo computazionale

n∑i=1

i +i−1∑j=1

j

+ n2 =

n∑i=1

i +n∑

i=1

i(i− 1)2

+ n2

=

n∑i=1

i +12

n∑i=1

i2 − 12

n∑i=1

i + n2

=12

n∑i=1

i +12

n∑i=1

i2 + n2

=n(n + 1)

4+

(n + 1)n(2n + 1)12

+ n2

=n3

6+

32

n2 +n3' n3

6




E possibile provare che l’algoritmo di Cholesky e stabile.

Ricordiamo che, peraltro, anche l’algoritmo di Gauss, senza pivoting,e stabile per le matrici simmetriche definite positive.




Molte routine automatiche calcolano R = LT anzicche L, ma eevidente che l’algoritmo e lo stesso (con un attento uso degli indici),data la simmetria della matrice di partenza A.

In tali casi dunque scriveremo

A = RTR.

Osserviamo infine che se si vuole risolvere il sistema

Ax = b

mediante la fattorizzazione di Cholesky bastera, una volta computatala fattorizzazione, e quindi calcolata L ( o rispettivamente R), risolverei seguenti due sistemi

Ly = b, LTx = y (RTy = b, Rx = y rispettivamente)




Le function Matlab

Chiudiamo questa carrellata dei principali metodi diretti per larisoluzione di un sistema lineare descrivendo le principali function diMatlab che implementano tali metodi.

Sia allora A una matrice non singolare di ordine n e b un vettorecolonna di ordine n.

risoluzione dei sistemi diagonali: se la matrice dei coefficienti Ae diagonale, A\b risolve il sistema mediante n operazioni didivisione;

risoluzione dei sistemi triangolari: se la matrice dei coefficienti Ae triagolare (superiore o inferiore) A\b risolve il sistemamediante algoritmo di sostituzione (all’indietro o in avanti aseconda struttura della matrice);




risoluzione di sistemi con matrice dei coefficienti qualsiasi: A\brisolve il sistema mediante metodo di eliminazione di Gauss (inrealta con la fattorizzazione LU) con pivoting e gli algoritmi disostituzione;

risoluzione di un sistema con matrice simmetrica e definitapositiva: A\b risolve il sistema mediante metodo di Cholesky egli algoritmi di sostituzione;

fattorizzazione LU: il comando [L,U,P] = lu(A) calcola ifattori L, U e la matrice di permutazione P tali che PA = LU;

fattorizzazione di Cholesky: se la matrice A e simmetrica edefinita positiva, il comando R = chol(A) calcola il fattoretriangolare superiore R tale che A = RTR;




calcolo del determinante di una matrice: det(A) calcola ildeterminante della matrice non singolare A, utilizzando lafattorizzazione LU;

condizionamento: il comando cond(A,p), con p = 1, 2, inf ,‘fro’, calcola il numero di condizionamento rispettivamente innorma 1, 2, infinito e Frobenius;

inversa di una matrice: inv(A), calcola l’inversa della matriceutilizzando la fattorizzazione LU.


Il metodo di Eliminazione di Gauss con la tecnica del pivoting

Documents

Transcript of Il metodo di Eliminazione di Gauss con la tecnica del pivoting