Metodi numerici per la risoluzione di Sistemi...

Metodi direttiMetodi Iterativi

Metodi numerici per la risoluzione di SistemiLineari

Stefano Berrone

Dipartimento di Matematicatel. 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


Problemi triangolariEliminazione GaussianaFattorizzazione LUPivoting

Problemi Triangolari - Sostituzione in avanti

Vogliamo risolvere il sistema lineare Ly = b con L triangolareinferiore (lower triangular) non singolare.

l11 0 0 ... 0l21 l22 0 ... 0l31 l32 l33 ... 0...

...ln1 ln2 ln3 ... lnn

y1

y2

y3...

yn

=

b1

b2

b3...

bn

Osservazione

Se L e non singolare, tutti gli elementi diagonali sono diversi da 0:

det(L) =n∏

i=1

lii




Esplicitando le equazioni:l11y1 = b1,l21y1 + l22y2 = b2,...ln1y1 + ln2y2 + ...+ lnnyn = bn.

=⇒ y1 = b1/l11

=⇒ y2 = (b2 − l21y1)/l22

=⇒ yi = (bi −i−1∑j=1

lijyj)/lii




Metodo di sostituzione in avanti (o forward substitution):

for i=1:ny(i)=b(i);for j=1:i-1

y(i)=y(i)-L(i,j)*y(j);endy(i)=y(i)/L(i,i);

end




Un’alternativa che sfrutta la capacita di MATLAB di operare susottomatrici e la seguente:

y(1) =b(1)/L(1,1);for i=2:ny(i)=(b(i)-L(i,1:i-1)*y(1:i-1))/L(i,i);end




Sostituzione all’indietro

Vogliamo risolvere il sistema lineare Ux = y con U triangolaresuperiore (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 nonsingolarita di U segue che uii 6= 0∀i




u11x1 + u12x2 + ... + u1nxn = y1,u22x2 + ... + u2nxn = y2,

...un−1,n−1xn−1 + un−1,nxn = yn−1,

unnxn = yn.

=⇒ xn = yn/unn.

=⇒ xn−1 = (yn−1 − un−1,nxn)/un−1,n−1.

=⇒ xi = (yi −n∑

j=i+1

uijxj)/uii .




Metodo di sostituzione all’indietro (o backward substitution):

for i=n:-1:1x(i)=y(i);for j=i+1:n

x(i)=x(i)-U(i,j)*x(j);endx(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




Eliminazione Gaussiana

La risoluzione di sistemi triangolari e quindi un problemaparticolarmente semplice da trattare.Idea: Nel caso di sistemi non triangolari, cercare di ricondursi aduno o piu sistemi equivalenti di forma triangolare.Dato un sistema lineare Ax = b le seguenti operazioni conduconoad 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 dellastessa e di un’altra equazione.




Metodo di eliminazione di Gauss:

consiste nell’applicare ripetutamente l’operazione 3 in modo da:

1 trasformare il sistema assegnato in un sistema a matricetriangolare superiore

2 e quindi nel risolvere il sistema equivalente cosı ottenuto conuna backward substitution.




Scriviamo il sistema Ax = b esplicitamente:

a11x1 + a12x2 + ... + a1nxn = b1,a21x1 + a22x2 + ... + a2nxn = b2,a31x1 + a32x2 + ... + a3nxn = b3,

...an1x1 + an2x2 + ... + annxn = bn.

Supponiamo a11 6= 0

r2 ← r2 + (−a21

a11) ∗ r1




a21x1 + a22x2 + ... + a2nxn = b2 +

−a21a11

( a11x1 + a12x2 + ... + a1nxn = b1 )

0x1 + a(2)22 x2 + ... + a

(2)2n xn = b

(2)2 .

a(2)22 = a22 − a12a21/a11, ...,a

(2)2n = a2n − a1na21/a11,

b(2)2 = b2 − b1a21/a11

Nella seconda equazione non compare piu la prima incognita.




Procediamo con tutte le ultime n − 1 equazioni:

ri ← ri + mi1 ∗ r1, i = 2, ..., n

mi1 = − ai1

a11, i = 2, ..., n.

=⇒

a(1)11 x1 + a

(1)12 x2 + ... + a

(1)1n xn = b

(1)1

0 + a(2)22 x2 + ... + a

(2)2n xn = b

(2)2

0 + a(2)32 x2 + ... + a

(2)3n xn = b

(2)3

...

0 + a(2)n2 x2 + ... + a

(2)nn xn = b

(2)n

a(1)ij = aij , b

(1)i = bi , i , j = 1, ..., n;

a(2)ij = a

(1)ij + mi1a

(1)1j , b

(2)i = b

(1)i + mi1b

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




Seconda colonna:

supponiamo a(2)22 6= 0. Eliminiamo x2 dalle ultime n − 2 equazioni

trasformate:

ri ← ri + r2 ∗mi2, mi2 = −a

(2)i2

a(2)22

i = 3, ..., n

=⇒

a(1)11 x1 + a

(1)12 x2 + a

(1)13 x3 + ... + a

(1)1n xn = b

(1)1

0 + a(2)22 x2 + a

(2)23 x3 + ... + a

(2)2n xn = b

(2)2

0 + 0 + a(3)33 x3 + ... + a

(3)3n xn = b

(3)3

...

0 + 0 + a(3)n3 x3 + ... + a

(3)nn xn = b

(3)n

a(3)ij = a

(2)ij + mi2a

(2)2j , b

(3)i = b

(2)i + mi2b

(2)2 , i , j = 3, ..., n




Dopo n − 1 passi:

a(1)11 x1 + a

(1)12 x2 + a

(1)13 x3 + ... + a

(1)1n xn = b

(1)1

0 + a(2)22 x2 + a

(2)23 x3 + ... + a

(2)2n xn = b

(2)2

0 + 0 + a(3)33 x3 + ... + a

(3)3n xn = b

(3)3

...

0 + 0 + 0 + ... + a(n)nn xn = b

(n)n

Osservazione

Il termine noto b ha subito esattamente tutte le trasformazionisubite dagli elementi della matrice A.




Definizione (pivot e moltiplicatori)

L’elemento a11 e i successivi elementi a(i)ii sono detti elementi

pivot. I coefficienti mij sono detti moltiplicatori.




Metodo di eliminazione di Gauss e fattorizzazione LU

In notazione matriciale, il k-esimo passo di eliminazione Gaussianapuo essere rappresentato come segue:

A(k+1) = MkA(k)

con

Mk =

1

1mk+1,k 1

.... . .

mn,k 1




Dopo n − 1 passi:

A(n) = Mn−1A(n−1) = Mn−1Mn−2A(n−2) = . . . = Mn−1 . . .M2M1︸︷︷︸M

A

b(n) = Mn−1 . . .M2M1b

eA(n)x = b(n)

con A(n) triangolare superiore

⇓Ax = b ⇔ MAx = Mb

Posto U = MA, si puo scrivere

Ux = Mb

da cui ottengo la soluzione con sostituzione all’indietro.Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari



Esempio

Dati

A =

4 4 82 8 71 3 6

b =

1297

calcolare la soluzione con eliminazione gaussiana.

M1 =

1 0 0−1

2 1 0−1

4 0 1

e

M1 ∗ (A|b) =

4 4 8 | 120 6 3 | 30 2 4 | 4

.




Esempio

M2 =

1 0 00 1 00 −1

3 1

e

M2M1 ∗ (A|b) =

4 4 8 | 120 6 3 | 30 0 3 | 3

= (U|b(3)).

Con sostituzione all’indietro:

x3 = 1

x2 = (3− 3 ∗ 1)/6 = 0

x1 = (12− 4 ∗ 0− 8 ∗ 1)/4 = 1




Com’e fatta M?- triangolare inferiore, ma

M 6=

1 0 0 . . . 0

m21 1 0 . . . 0m31 m32 1 . . . 0

......

......

...mn1 mn2 mn3 . . . 1

Poniamo L = M−1

U = MA =⇒ LU = LMA = A

Quindi il sistema di partenza equivale a

LUx = b.

Com’e fatta L? Suo ruolo?Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari



Struttura di L

L triangolare inferiore.Il calcolo e meno complicato e costoso di quanto sembri a primavista: non occorre costruire M e poi calcolarne l’inversa.

L =

1 0 0 . . . 0−m21 1 0 . . . 0−m31 −m32 1 . . . 0

.... . .

−mn1 −mn2 . . . −mn,n−1 1

Osservazione

L puo quindi essere costruita durante l’eliminazione gaussianasemplicemente conservando gli elementi mij e cambiandogli ilsegno.




Fattorizzazione LU

Abbiamo fattorizzato A = LU con

L triangolare inferiore

U triangolare superiore

Utilizzo per risolvere il sistema lineare:

Ax = b

LUx = b

L Ux︸︷︷︸y

= b

⇓

Ly = b, Ux = y




Il problema di partenza puo essere affrontato in due fasi, risolvendodue 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 e del tutto equivalente adeffettuare sul vettore b le trasformazioni dell’eliminazionegaussiana. Quindi il vettore y calcolato come soluzione del sistemalineare Ly = b coincide con b(n)!




Meglio fattorizzazione LU + 2 sistemi lineari triangolari, oeliminazione gaussiana senza costruire esplicitamente L?Con l’eliminazione Gaussiana, b subisce le stesse trasformazioni diA. Su un solo sistema sono equivalenti. Se devo risolvere unsecondo sistema lineare

Ax = c

con eliminazione gaussiana devo rieseguire tutto ilprocedimento per effettuare le trasformazioni su c e ottenerec(n);

con fattorizzazione LU gia calcolata devo solo risolvere isistemi Ly = c e Ux = y .

E veramente piu conveniente fare cosı piuttosto che rieffettuareex-novo l’eliminazione gaussiana sul sistema Ax = c?




Costo computazionale

Per costo computazionale di un algoritmo intendiamo l’ordine digrandezza del numero di moltiplicazioni o divisioni effettuate,calcolato in funzione della dimensione n del problema (in questocaso l’ordine della matrice).

Costo computazionale per calcolo fattorizzazione LU(coincidente con quello della sola eliminazione gaussiana):

O(

n3

3

)moltiplicazioni.

Il costo computazionale soluzione di un sistema triangolare(superiore o inferiore indifferentemente):

O(

n2

2

)moltiplicazioni.




Risolvendo un sistema lineare con fattorizzazione LU la maggiorparte del costo risiede nella fattorizzazione e non nella soluzionedei sistemi triangolari ottenuti: il costo totale e infatti di

O(n3

3) +O(n2) ' O(

n3

3)

operazioni, essendo (per n grande) n2 trascurabile rispetto a n3.

Osservazione

I metodi di fattorizzazione sono quindi particolarmente convenientiquando si debbano risolvere un certo numero di sistemi lineari incui la matrice dei coefficienti e sempre la stessa e cambia solo ilvettore dei termini noti.




Esempio

Data la matrice

A =

4 4 82 8 71 3 6

calcolarne la fattorizzazione LU.Basta recuperare i moltiplicatori:

U =

4 4 80 6 30 0 2

L =

1 0 012 1 014

13 1

Osservazione: det(A) = det(L) · det(U) = 1 ·

∏ni=1 uii = 48.




Osservazione

L puo essere memorizzata nella stessa matrice A via via che inquesta vengono azzerati gli elementi sotto la diagonale.

4 4 82 8 71 3 6

−→4 4 8

12 6 314 2 4

−→4 4 8

12 6 314

13 2




L’eliminazione di Gauss / fattorizzazione LU puo essere sempreportata a termine con successo?

Definizione

Sia A ∈ Rn×n; definiamo Ak , con k = 1, ..., n, sottomatriceprincipale di testa di ordine k la matrice ottenuta intersecando leprime k righe e k colonne di A.

Teorema

Sia A una matrice di ordine n. Se Ak e non singolare perk = 1, ..., n − 1 allora esiste ed e unica la fattorizzazione LU di A.




Osservazione

Le matrici a diagonale dominante per righe o per colonne e lematrici simmetriche definite positive soddisfano le ipotesi delteorema precedente.

Osservazione

Il teorema precedente non esclude che la matrice A sia singolare.Se anche A e singolare, ma det(Ak) 6= 0, k = 1, ..., n − 1, esistecomunque la fattorizzazione LU, con la matrice U che avra unelemento uii nullo.




La Strategia del Pivoting

Esempio

Che sucede se voglio applicare eliminazione gaussiana a

A =

(0 11 0

)?

Occorre ricorrere alle matrici di permutazione.




Definizione

Sia Pij la matrice ottenuta dalla matrice identita permutando lai-esima e la j-esima riga:

Pij =

i j

1 0 0 . . . 0 . . . 0 . . . 00 1 0 . . . 0 . . . 0 . . . 00 0 1 . . . 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



Data una matrice A

la moltiplicazione a sinistra di A per Pij (PijA) produce loscambio della i-esima e j-esima riga di A;

la moltiplicazione a destra di A per Pij (APij) produce loscambio della i-esima e j-esima colonna di A.

Se moltiplichiamo piu matrici di permutazione semplice fra loro,otteniamo una matrice (che indicheremo con P e che chiameremomatrice di permutazione) il cui effetto su una matrice A saraquello di effettuare tutti gli scambi di righe (colonne) associate allematrici di permutazione semplice.




Teorema

Sia A una matrice di ordine n. Allora esiste una matrice dipermutazione P per cui si puo ottenere la fattorizzazione LU di A,cioe:

PA = LU.

Si osservi che, data A, vi possono essere diverse matrici dipermutazione P tali che PA ammette fattorizzazione LU.

Come individuare una matrice di permutazione?

⇓

Strategia di pivoting parziale.




Pivoting Parziale

Se durante l’eliminazione di Gauss al passo k-esimo risulta

a(k)kk = 0, scambiamo due equazioni in modo che il nuovo

elemento pivot a(k)kk = 0 sia diverso da zero.

La scelta e sempre tra le equazioni che seguono la k-esima.




Esempio

Consideriamo il seguente sistema lineare:x1 + x2 + x3 = 1x1 + x2 + 2x3 = 2x1 + 2x2 + 2x3 = 1

Eliminazione gaussiana: m21 = −1 e m31 = −1 =⇒x1 + x2 + x3 = 10 + 0 + x3 = 10 + x2 + x3 = 0

a(2)22 = 0!




Esempio (segue)

;

x1 + x2 + x3 = 10 + x2 + x3 = 00 + 0 + x3 = 1

Questo e gia in forma triangolare e si puo ottenere per sostituzioneall’indietro la soluzione:

x3 = 1x2 = −x3 = −1x1 = 1− x2 − x3 = 1

Abbiamo det(A) = − det(U) = −1; in generale,det(A) = (−1)s · det(U), dove s e il numero di permutazioni diriga effettuate.




Fattorizzazione PA = LU

Cosa accade alla fattorizzazione LU con il pivoting?

A(k+1) = MkPkA(k)

n − 1 passi:

A(n) = Mn−1Pn−1A(n−1) = ... = Mn−1Pn−1Mn−2Pn−2 . . .M2P2M1P1A

U = Mn−1Pn−1Mn−2Pn−2Mn−3 . . .M2P2M1P1A

= Mn−1Mn−2Pn−1Pn−2Mn−3 . . .M2P2M1P1A

= Mn−1Mn−2 . . . M2M1︸︷︷︸M

Pn−1Pn−2 . . .P2P1︸︷︷︸P

A

= MPA




U = MPA

M−1U = 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




In generale la strategia di permutare equazioni/righe viene usataanche se si incontrano elementi pivot diversi da zero peraumentare la stabilita dell’eliminazione gaussiana.Si puo infatti dimostrare che per aumentare la stabilita e opportunoevitare moltiplicatori troppo grandi, poiche amplificano gli erroridi arrotondamento, con conseguente instabilita dell’algoritmo.Per prevenire la crescita dei moltiplicatori e necessario prevenireelementi pivot troppo piccoli.Per raggiungere questo obiettivo possiamo applicare la strategiadel pivoting parziale ammettendo scambi di righe anche se al

passo k l’elemento pivot a(k)kk non e nullo.

Piu precisamente, al passo k si determina l’elemento di modulo

massimo fra tutti gli elementi a(k)rk , con r ≥ k; una volta

individuato l’indice r corrispondente all’elemento di modulomassimo si scambiano fra loro la riga k-esima e r -esima.




In questo modo ci si assicura che ad ogni passo si sta usandol’elemento pivot piu grande possibile. Questa strategiaincrementa in modo significativo la stabilita dell’algoritmo.

Osservazione

Se la matrice A e a diagonale dominante per colonne e possibiledimostrare che la strategia del pivoting parziale e del tuttosuperflua in quanto non provoca alcun scambio diequazioni.Inoltre, se la matrice A e simmetrica e definita positiva, ilmetodo di Gauss senza pivoting risulta numericamente stabile,quindi non c’e bisogno di effettuare pivoting. Per tale classe dimatrici esiste una diversa fattorizzazione molto piu conveniente dausare (la fattorizzazione di Cholesky).




Come costruire L

P = Pn−1 . . .P1

U = A(n)

L =?

E facilmente ottenibile con l’eliminazione gaussiana tenendo contodelle permutazione anche nell’immagazzinamento dei moltiplicatorisotto la diagonale.




Esempio

A =

1 4 82 0 74 2 6

P1 = P13:

P13A =

4 2 62 0 71 4 8

A(2) =

4 2 612 −1 414

72

132




Esempio (segue...)

P2 = P23:

P23A(2) =

4 2 614

72

132

12 −1 4

A(3) =

4 2 614

72

132

12 −2

7417

L =

1 0 014 1 012 −2

7 1

U =

4 2 60 7

2132

0 0 417

P = P23P13 =

0 0 11 0 00 1 0




Fattorizzazione di Cholesky

Se A e simmetrica e definita positiva, essendo l’eliminazionegaussiana stabile senza pivoting, si puo sfruttare la simmetria di Aper ridurre il costo computazionale del calcolo della fattorizzazioneusando la fattorizzazione di Cholesky anziche fattorizzazione LU.

Teorema (Fattorizzazione di Cholesky)

Sia A una matrice di ordine n simmetrica e definita positiva. Alloraesiste un’unica matrice triangolare superiore R con elementidiagonali positivi tale che

A = RT R.

Costo computazionale per fattorizzazione di Cholesky:

O(n3

6)




Comandi MATLAB

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

2 R=chol(A) Attenzione! chol assume che A sia simmetricadefinita positiva e usa solo la parte triangolare superiore!

3 \

Altre fattorizzazioni

1 Fattorizzazione QR (qr(A)): Q e ortogonale, R e triangolaresuperiore. Molto stabile ma anche piu costosa dellafattorizzazione PA = LU. Utile per altre applicazioni.

Varie

1 norm(x), norm(A): calcolo di norme

2 cond(A): numero di condizionamento



GeneralitaMetodo di JacobiMetodo di Gauss-SeidelTest di arrestoCondizioni 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 dipassi.Per sistemi con matrici dense i metodi diretti sono generalmentepreferibili.Per sistemi con matrici sparse e di ordine elevato, i metodi direttidiventano impraticabili a causa del fill in ⇒ non puo esseresfruttata la sparsita della matrice.




I metodi iterativi lasciano inalterata la struttura di A ecostruiscono una successione infinita di vettori {x (k)}k≥0, chesotto opportune ipotesi converge alla soluzione x .Per i metodi iterativi dovremo discutere:

1 condizioni di applicabilita

2 costo computazionale (di ogni iterazione)

3 problemi connessi alla convergenza: convergenza si/no evelocita di convergenza.




Definizione

Una successione di vettori {x (k)}k≥0 di Rn si dice convergente adun vettore x ∈ Rn se esiste una norma per cui

limk→∞

‖x (k) − x‖ = 0,

ed in tal caso si ponelim

k→∞x (k) = x .

Osservazione

Grazie all’equivalenza delle norme su Rn tale definizione risultaindipendente dalla norma scelta in Rn; se una successione di vettoriconverge in una norma di Rn, allora converge in tutte le norme diRn.




Osservazione

La condizione di convergenza in una qualunque norma di Rn sitraduce in una condizione di convergenza per componenti:

limk→∞

x (k) = x ⇐⇒ limk→∞

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

limk→∞

Ak = 0 ∈ Rn×n ⇐⇒ ρ(A) < 1

Ricordiamo che per ogni norma di matrice compatibile con unanorma di vettore si ha ρ(A) ≤ ‖A‖.




Generalita sui metodi iterativi

Idea base: data una stima iniziale x (0) della soluzione del sistemalineare Ax = b, con A ∈ Rn×n, det(A) 6= 0, b ∈ Rn, si costruisceuna successione di vettori {x (k)}k≥0 che converga (si spera) ad x ,risolvendo ad ogni iterazione dei sistemi lineari molto piu semplicidi 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.




Supponiamo di essere all’iterazione k + 1, di avere quindi adisposizione 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−1Nx (k) + M−1b

che puo essere scritto in forma compatta come

x (k+1) = Bx (k) + c ,

avendo posto B = −M−1N e c = M−1b. La matrice

B = −M−1N

e la matrice d’iterazione ed individua il particolare metodo; efondamentale nello stabilire le condizioni e la velocita diconvergenza.




Proprieta (consistenza)

Se x (k) → x , allora x soddisfa l’equazione

x = Bx + c

e quindi il limite della successione x e proprio la soluzione x diMx = −Nx + b (cioe 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 inizialex (0) la successione {x (k)}k≥0 e convergente.




Studio della convergenza dei metodi iterativi

Sotto quali condizioni la successione {x (k)}k≥0 generata da unmetodo iterativo converge?Se la successione converge ad un vettore x , questo e senz’altrosoluzione del sistema lineare. Definiamo l’errore e(k) commesso alpasso k come la differenza fra l’approssimazione corrente e lasoluzione:

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)) = B2e(k−2) = ... = Bke(1) = Bke(0).Stefano Berrone Metodi numerici per la risoluzione di Sistemi Lineari



Poichee(k) = Bke(0)

abbiamo

limk→∞

x (k) = x ⇔ limk→∞

e(k) = 0⇔ limk→∞

Bke(0) = 0⇔ limk→∞

Bk = 0

Quindi c’e convergenza se e solo se ρ(B) < 1.




Teorema

Il metodo iterativo definito dalla matrice di iterazione B convergeper ogni x (0) ∈ Rn se e solo se

ρ(B) < 1.

La condizione necessaria e sufficiente espressa dal teorema per laconvergenza non e di facile verifica. Una condizione sufficiente eespressa dal seguente Teorema.

Teorema

Sia ‖ · ‖ una norma di matrice naturale. Se ‖B‖ < 1 il metodoiterativo converge.Dimostrazione: E sufficiente ricordare che per ogni norma naturale

ρ(B) ≤ ‖B‖




Il raggio spettrale indica non solo se il metodo iterativoconverge ma anche quanto velocemente converge.

Intuitivamente, la velocita di convergenza e una misura diquanto velocemente l’errore tende a zero.

La velocita di convergenza di un metodo iterativo e tantomaggiore quanto piu il raggio spettrale della matrice diiterazione e vicino a 0; viceversa e tanto piu lenta quanto piue vicino a 1.




Metodo di Jacobi

Detto anche metodo delle sostituzioni simultanee.Consideriamo il seguente splitting di A:

A = E + D + F , A =

FD

E

dove

1 E e la parte di A sotto la diagonale,

2 D e la diagonale di A,

3 F e la parte di A sopra la diagonale

Prendiamo M = D, N = E + F . Matrice d’iterazione BJ :

BJ = −D−1(E + F )




Data l’iterata corrente x (k) ∈ Rn, la x (k+1) ∈ Rn viene calcolatatramite il seguente algoritmo:

for i = 1 : nx

(k+1)i =

(bi −

∑i−1j=1 aijx

(k)j −

∑nj=i+1 aijx

(k)j

)/aii

endIn forma compatta, puo essere scritto

x (k+1) = D−1(

b − Ex (k) − Fx (k)).

Osservazione

Se aii = 0 per qualche i , D non e invertibile. Si possono peropreliminarmente riordinare le equazioni del sistema lineare in mododa portare in posizione ii un elemento non nullo. Cio e semprepossibile se A e non singolare.




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)−1F




Data l’iterata corrente x (k) ∈ Rn, la x (k+1) ∈ Rn viene calcolatatramite il seguente algoritmo:

for i = 1 : nx

(k+1)i =

(bi −

∑i−1j=1 aijx

(k+1)j −

∑nj=i+1 aijx

(k)j

)/aii

endIn forma compatta, puo essere scritto

x (k+1) = D−1(

b − Ex (k+1) − Fx (k)).




Test di arresto

Un metodo iterativo fornisce una successione infinita di iterate x (k).Questo diventera un algoritmo solo se stabiliremo in modo precisoquando interrompere il procedimento di generazione di iterate.Dobbiamo ovvero stabilire un criterio di arresto delle iterazioni.

1 si arresta il processo iterativo quando risulta∥∥ x (k+1) − x (k)∥∥∥∥ x (k+1)

∥∥ ≤ tol,

dove tol e una tolleranza relativa fissata sulla base diopportuni criteri di utilizzo della soluzione.

2 oppure se ∥∥∥ x (k+1) − x (k)∥∥∥ ≤ tol,

dove tol e questa volta una tolleranza assoluta.




Una strada alternativa per definire un test d’arresto consiste nelcontrollo del residuo dell’equazione. Sia x (k+1) la soluzioneottenuta al passo k + 1, il residuo corrispondente e definito come

r (k+1) = b − Ax (k+1).

Test di arresto relativo sul residuo:∥∥ b − Ax (k+1)∥∥

‖ b ‖≤ tol.




Osservazione

Residuo ”piccolo” non garantisce errore ”piccolo”.∥∥∥ e(k)∥∥∥ =

∥∥∥ x (k) − x∥∥∥ =

∥∥∥ x (k) − A−1b∥∥∥ =

∥∥∥A−1(

Ax (k) − b)∥∥∥

=∥∥∥−A−1r (k)

∥∥∥ ≤ ∥∥A−1∥∥∥∥∥ r (k)

∥∥∥ .Ricordando 1

‖ x ‖ ≤‖A ‖‖ b ‖ (v. conti su condizionamento) abbiamo∥∥ e(k)

∥∥‖ x ‖

≤∥∥A−1

∥∥∥∥∥ r (k)∥∥∥ ‖A ‖‖ b ‖

= K (A)

∥∥ r (k)∥∥

‖ b ‖

Quindi∥∥ r (k)∥∥

‖ b ‖≤ tol =⇒

∥∥ e(k)∥∥

‖ x ‖≤ K (A)

∥∥ r (k)∥∥

‖ b ‖≤ K (A)tol

Osservazione

In ogni caso, qualsiasi test di arresto si usi e sempre indispensabileprevedere un numero massimo di iterazioni.




Teorema (Condizioni di convergenza per il metodo di Jacobi)

Per matrici a predominanza diagonale stretta (per righe o percolonne 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 percolonne) 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-Seidelconverge.




Confronto delle proprieta di convergenza del metodo diJacobi e di Gauss-Sidel

In generale la convergenza di uno dei metodi non implica laconvergenza dell’altro e viceversa.Tuttavia quando entrambi convergono la velocita di convergenzadel metodo di Gauss-Seidel e generalmente superiore.

Teorema (Teorema di Stein-Rosenberg)

Sia A ∈ Rn×n con aij ≤ 0, ∀i 6= j e aii > 0, allora si verifica uno euno solo dei seguenti risultati:

1 0 < ρ(BGS) < ρ(BJ) < 1,

2 1 < ρ(BJ) < ρ(BGS),

3 ρ(BGS) = ρ(BJ) = 0,

4 ρ(BGS) = ρ(BJ) = 1.




Teorema

Sia A ∈ Rn×n una matrice tridiagonale con elementi diagonali nonnulli, allora

ρ(BGS) = ρ2(BJ),

cioe il metodo di Gauss-Seidel ed il metodo di Jacobi convergono odivergono simultaneamente e il tasso asintotico di convergenza delmetodo di Gauss-Seidel e doppio di quello del metodo di Jacobi.


Metodi numerici per la risoluzione di Sistemi...

Documents

Transcript of Metodi numerici per la risoluzione di Sistemi...