SISTEMI LINEARI QUADRATI: METODI ITERATIVI

SISTEMI LINEARI QUADRATI:METODI ITERATIVI

CALCOLO NUMERICO e PROGRAMMAZIONE

SISTEMI LINEARI QUADRATI:METODI ITERATIVI – p.1/54

RICHIAMI di ALGEBRA LINEAREDEFINIZIONI

A ∈ Rn×n simmetrica se A = AT ;

A ∈ Cn×n Hermitiana se A = A∗, dove A∗ = (aj,i) i=1,...,nj=1,...,n

;

A ∈ Rn×n ortogonale se AT A = I;

A ∈ Cn×n (simmetrica) unitaria se A∗A = I;

A ∈ Rn×n (simmetrica) definita positiva se xT Ax > 0, ∀x 6= 0;

A ∈ Cn×n (hermitiana) definita positiva se x∗Ax > 0, ∀x 6= 0;

A diagonale dominante in senso forte[debole] se

|ai,i| >

n∑

j=1

j 6=i

|ai,j |,

|ai,i| ≥

n∑

j=1

j 6=i

|ai,j |

i = 1, . . . , n

A diagonale dominante se

|ai,i| ≥n

∑

j=1

j 6=i

|ai,j |, i = 1, . . . , n, ed ∃k : |ak,k| >

n∑

j=1

j 6=k

|ak,j |


RICHIAMI di ALGEBRA LINEARETEOREMASe A ∈ Rn×n è simmetrica, sono equivalenti:

1. A definita positiva;

2. gli autovalori di A sono reali e positivi;

3. i determinanti dei minori principali di A verificano: det(Ak) > 0, k = 1, . . . , n

(criterio di Sylvester).

L’algebra lineare numerica richiede spazi lineari normati

ESEMPIO: Si consideri il sistema lineare

8 1 1

1 5 −1

1 −1 5

x1

x2

x3

=

26

7

7

la cui soluzione e x = [3, 1, 1]T . Si supponga di aver ottenuto le seguenti approssimazioni della

soluzione: v = [2.9997675, 1.0002365, 1.0003005]T

z = [3.0002725, 0.9996165, 0.9995525]T

Come possiamo stabilire quale delle due grandezze e piu vicina alla soluzione del sistema?


Norma di vettoreDefinizione: Sia X spazio lineare e f : X → R funzionale tale che

1. f(x) ≥ 0, ∀x ∈ X, f(x) = 0 ⇔ x = 0;

2. f(x + y) ≤ f(x) + f(y), ∀x, y ∈ X;

3. f(αx) = |α|f(x), ∀α ∈ R(C), ∀x ∈ X;

X è detto spazio lineare normato (sln) e f è detta norma.In particolare per X = Rn si usano le norme p o Hölderiane:

‖x‖p := (|x1|p + |x2|p + · · · + |xn|p)1

p

‖x‖1 := |x1| + |x2| + · · · + |xn|

‖x‖2 :=√

x21 + x2

2 + · · · + x2n

‖x‖∞ := max{|x1|, |x2|, . . . , |xn|}


Norma di vettore

Definizione: Siano f ed g due norme di Rn. Se esistono c, C ∈ R+ t.c.

cg(x) ≤ f(x) ≤ Cg(x), ∀x ∈ Rn

f ed g si dicono norme equivalenti.

‖ · ‖1, ‖ · ‖2, ‖ · ‖∞ sono norme equivalenti

‖x‖∞ ≤ ‖x‖1 ≤ n‖x‖∞‖x‖∞ ≤ ‖x‖2 ≤

√n‖x‖∞

‖x‖1/√

n ≤ ‖x‖2 ≤ ‖x‖1

ma con diversa geometria: S = {x ∈ Rn : ‖x‖ ≤ 1}

|| ||1

|| ||2

|| ||∞


Norma di matrici in Rn×n

Definizione: Sia ‖ · ‖Rn una norma vettoriale. ‖ · ‖Rn×n si dice compatibile o

consistente con ‖ · ‖Rn se ∀A, B ∈ Rn×n e ∀x ∈ Rn si ha

‖Ax‖Rn ≤ ‖A‖Rn×n‖x‖Rn

‖AB‖Rn×n ≤ ‖A‖

Rn×n‖B‖Rn×n

Definizione: Sia ‖ · ‖Rn×n compatibile, si dice indotta da ‖ · ‖Rn se ∀A ∈ Rn×n

∃x ∈ Rn t.c.

‖Ax‖Rn = ‖A‖Rn×n‖x‖Rn

Definizione: Sia data ‖ · ‖Rn norma vettoriale, il funzionale ‖ · ‖ : Rn×n → R t.c.

‖A‖ := supx6=0

‖Ax‖Rn

‖x‖Rn

= max‖x‖Rn=1

‖Ax‖Rn

è una norma, norma naturale

Ogni norma naturale è indotta:sia y : ‖Ay‖ = max‖x‖=1 ‖Ax‖ ⇒ ‖Ay‖ = ‖A‖ = ‖A‖‖y‖


Principali norme indotte naturali

‖ ‖1 induce ‖A‖1 := max1≤j≤n

n∑

i=1

|ai,j |

‖ ‖∞ induce ‖A‖∞ := max1≤i≤n

n∑

j=1

|ai,j |

‖ ‖2 induce ‖A‖2 :=√

ρ(AT A)

dove data M ∈ Rp×p, ρ(M) = max1≤k≤p |λk(M)| è il suo raggio spettrale


Proprietà delle norme naturali

Sia M ∈ Rp×p

ρ(M) = inf{‖M‖ : ∀‖ · ‖ norma naturale}

ovvero ∀ε > 0, ∃‖ · ‖ norma naturale t.c ρ(M) ≤ ‖M‖ ≤ ρ(M) + ε, ne segue

∀M ∈ Rp×p, ∀‖ · ‖ norma naturale,

ρ(M) ≤ ‖M‖

Una matrice M ∈ Rp×p si dice convergente se limk→∞ Mk = 0

Teorema

M convergente ⇔ ρ(M) < 1

Invarianza della norma rispetto a trasformazioni ortogonali:

x ∈ Rn, U ∈ Rn×n ortogonale: ‖x‖2 = ‖Ux‖2

A ∈ Rm×n, U ∈ Rm×m, V ∈ Rn×n ortogonali: ‖A‖2 = ‖UAV ‖2


SISTEMI LINEARI QUADRATI

Sia A ∈ Rn×n, b ∈ Rn, se det A 6= 0 allora ∃!x ∈ Rn t.c.

Ax = b.

Molti modelli matematici significativi sono di tipo lineare

Sistemi lineari nascono da contesti diversi, è importante disporre di un ampiavarietà di algoritmi in modo da scegliere il più adatto al problema specifico(stabilità, occupazione di memoria, velocità)

I metodi per i sistemi lineari si dividono in due gruppi:

Metodi diretti: si basano sull’idea di trasformare il sistema attraverso unnumero finito di operazioni in un sistema equivalente di cui siaesplicitamente calcolabile la soluzione. In assenza di errori diarrotondamento forniscono la soluzione esatta.

Metodi iterativi: la soluzione è ottenuta come limite di una successione.Permettono di sfruttare la sparsità della matrice dei coeffcienti in quanto alcontrario dei metodi diretti, tale matrice non viene modificata


CondizionamentoDato Ax = b si vuol studiare come varia la soluzione del sistema al variare deidati(A, b).

Si studia il caso in cui solo il termine noto viene modificato. Si esprime laperturbazione sui dati con δb ∈ Rn e siano x e x + δx rispettivamente lesoluzioni di

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

sottraendo i due sistemi

Aδx = δb ⇒ δx = A−1δb

da cui rispetto ad una norma naturale:

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

Poichè ‖b‖ = ‖Ax‖ ≤ ‖A‖ ‖x‖ ne segue 1/‖x‖ ≤ ‖A‖/‖b‖ e quindi

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

‖b‖


CondizionamentoConsideriamo il caso in cui la matrice A venga perturbata, si esprimeεδA ∈ Rn×n la perturbazione di A e εδb ∈ Rn quella di b (ε ∈ R positivo esufficientemente piccolo). Si risolve il sistema

(A + εδA)x(ε) = b + εδb; (1)

con x(0) = x. Se det A 6= 0 allora per ε “piccolo” det(A + εδA) 6= 0; posso ricavare

x(ε) = (A + εδA)−1(b + εδb)

funzione di ε derivabile in un intorno dell’origine.Derivando (1) rispetto a ε si ha

δAx(ε) + (A + εδA)x(ε) = δb;

per ε = 0 : δAx(0) + Ax(0) = δb da cui x(0) = A−1(δb − δAx(0)).Sviluppando in serie x(ε) = x(0) + εx(0) + O(ε2) ⇒ x(ε) − x(0) ' εx(0) da cui

‖x(ε) − x(0)‖‖x(0)‖ ' ‖εx(0)‖

‖x(0)‖ = ε‖A−1(δb − δAx(0))‖

‖x(0)‖


Condizionamento

‖x(ε) − x(0)‖‖x(0)‖ ≤ ε

‖A−1‖‖δb − δAx(0)‖‖x(0)‖

≤ ε‖A−1‖( ‖δb‖‖x(0)‖ + ‖δA‖

)

= ε‖A−1‖‖A‖( ‖δb‖‖A‖‖x‖ +

‖δA‖‖A‖

)

≤ ‖A−1‖‖A‖(‖εδb‖

‖b‖ +‖εδA‖‖A‖

)

≤ κ(A)

(

ε‖δb‖‖b‖ + ε

‖δA‖‖A‖

)

κ(A) := ‖A−1‖‖A‖: numero di condizionamento (‖ · ‖ qualsiasi normanaturale)

‖x(ε) − x‖‖x‖ ≤ κ(A)

(

ε‖δb‖‖b‖ + ε

‖δA‖‖A‖

)

Il rapporto tra l’errore relativo sui risultati e l’errore relativo sui dati risultaminore o uguale di κ(A)


Condizionamento: proprietà

κ(A) ≥ 1: infatti κ(A) = ‖A−1‖‖A‖ ≥ ‖A−1A‖ = ‖I‖ = 1

Non ci sono relazioni tra ordine del sistema e numero di condizionamento:

A1 =

−1 1

−1 1 + 10−5

, A−11 = 105

−1 1

−1 1

κ∞(A1) ' 4 · 105: matrice piccola e malcondizionata

A2 =

1 −1 −1 . . . −1

0 1 −1 . . . −1

0 0 1 . . . −1

.... . .

...

0 . . . 1

, A−12 =

1 1 2 4 . . . 2n−2

0 1 1 2 . . . 2n−3

0 0 1 . . . 2n−4

.... . .

...

0 . . . 1

κ∞(A2) = ‖A2‖‖A−12 ‖ = n · 2n−1: matrice grande e malcondizionata


Condizionamento: proprietà

Non ci sono relazioni tra determinante e numero di condizionamento:det(A2) = 1 ≫ 0, ma κ∞(A2) = n2n−1 malcondizionata

A3 =

110

0 . . . 0

0 110

. . . 0

.... . .

...

0 . . . 110

, A−13 =

10 0 . . . 0

0 10 . . . 0

.... . .

...

0 . . . 10

,

det(A3) = 10−n, κ∞(A3) = 1: matrice quasi singolare e bencondizionata

Non è semplice calcolare κ(A): occorre calcolare A−1 che equivale a risolveren sistemi lineari ⇒ è importante stimare κ(A).


Dato Ax = b, come calcolare x?

Il metodo di Cramer è inefficiente.

Esempio: n = 30 (sistema medio-piccolo)Occorrono n + 1 determinanti di ordine n,per ciascuno occorre sommare n! termini.In ognuno ci sono n − 1 prodotti, complessivamente occorrono:(n + 1)n!(n − 1) prodotti ovvero per n = 30 circa 1034 prodotti e supponendoun tempo di 10−6 secondi per prodotto, si ottiene un tempo di circa1028sec ' 3 · 1020 anni

In qualche caso il calcolo è semplice: sistemi con matrici triangolari inferiori osuperiori.

A =

2 0 0

2 4 0

3 −2 3

:

2x = 4

2x + 4y = 8

3x − 2y + 3z = 10

x = 2

y = 8−2·24

3x − 2y + 3z = 10

x = 2

y = 1

z = 10−3·2+2·13

= 2


Sistemi triangolari

Ly = b con L ∈ Rn×n triangolare inferiore non singolare

`11y1 = b1

`21y1 + `22y2 = b2

· · ·ì1y1 + · · · + ìiyi = bi

· · ·`n1y1 + · · · · · · + `nnyn = bn

⇒

y1 = b1/`11

y2 = (b2 − `21y1)/`22

· · ·yi = [bi − (ì1y1 + · · · + ì,i−1yi−1)]/ìi

· · ·yn = (bn − ∑n−1

j=1 `njyj)/`nn

Algoritmo di Sostituzione in avanti

1. Per i = 1, 2, . . . , n

1. yi := bi

2. Per j = 1, 2, . . . , i − 1

yi := yi − ìjyj

3. yi := yi/ìi


Sistemi triangolari

Rx = y con R ∈ Rn×n triangolare superiore non singolare

r11x1 + r12x2 · · · + r1nxn = y1

r22x2 + · · · + r2nxn = y2

...

rnnxn = yn

Algoritmo di Sostituzione all’indietro

1. Per i = n, n − 1, . . . , 1

1. xi := yi

2. Per j = i + 1, i + 2, . . . , n

xi := xi − rijxj

3. xi := xi/rii

costo computazionale: si calcola il numero di prodotti/divisioni eseguite infunzione di n e si assume come costo il termine con la potenza più alta:

n∑

i=2

(i − 1) =

n−1∑

i=1

i =n(n − 1)

2prodotti

n divisioni ⇒ n2

2operazioni


Metodi diretti: fattorizzazione LRIdea Base: trasformare il problema complesso in uno equivalente più semplice.

Se esistono L ∈ Rn×n triangolare inferiore e R ∈ Rn×n triangolare superioretali che A = LR (fattorizzazione LR) allora:

AX = b ⇔ LRx = b ⇔

Ly = b

Rx = y

Esiste la fattorizzazione LR e se esiste è unica?

`11 0 · · · 0

`21 `22. . .

......

. . . 0

`n1 `n2 · · · `nn

·

r11 r12 · · · r1n

0 r22 · · · r2n

.... . .

. . ....

0 · · · 0 rnn

=

a11 · · · a1n

......

......

an1 · · · ann

n2 equazioni non lineari in n2 + n incognite ⇒ si fissano n incognite

r11 = r22 = · · · = rnn = 1: fatt. alla Crout`11 = `22 = · · · = `nn = 1: fatt. alla Dolittle ( associata al metodo di

Gauss)SISTEMI LINEARI QUADRATI:METODI ITERATIVI – p.18/54

Metodi iterativiNei metodi diretti la presenza di eventuali elementi nulli nella matrice non puòessere sfruttata ai fini di ridurre in il costo computazionale sia l’occupazione dimemoria ( aspetti significativi per sistemi di grandi dimensioni)

infatti la trasformazione di A può introdurre un numero diverso da zero laddoveprima c’era uno zero (fill in)

I metodi iterativi sono utili per la risoluzione di sistemi lineari di grandidimensioni con matrici A sparse ( il numero degli elementi non nulli èdell’ordine di n)

Dato il sistema lineare

Ax = b, A ∈ Rn×n, b ∈ R

n

con soluzione x∗ ∈ Rn, si vuol costruire partendo, da un x(0) ∈ Rn assegnato, unasuccessione {x(k)} tale che

(1) x(k) → x∗

(2) x(k) facile e non dispendioso da costruire


Metodi iterativi

ESEMPIO: Assegnato x(0) ∈ Rn

METODO di JACOBI (o spostamenti simultanei)

x(k+1)i =

1

aii

bi −n

∑

j=1,j 6=i

aijx(k)j

aii 6= 0

METODO di GAUSS-SEIDEL (o spostamenti successivi)

x(k+1)i =

1

aii

bi −i−1∑

j=1,j 6=i

aijx(k+1)j −

n∑

j=i+1,j 6=i

aijx(k)j

aii 6= 0

si può notare che nella prima sommatoria si usano le componenti “nuove” delvettore


Metodi iterativiPer entrambi si può dare una formulazione matriciale: scomponiamo

A = L + D + U

L =

0 0 · · · 0

a21 0. . .

...

a31 a32 0. . .

......

. . . 0

an1 an2 · · · an,n−1 0

D =

a11 0 · · · 0

0 a22

......

...

0 · · · ann

U =

0 a12 · · · a1n

0 0 a23 · · · a2n

.... . .

. . ....

0 · · · 0 an−1,n

0 · · · 0 0

Jacobi: aiix(k+1)i = −∑i−1

j=1,j 6=iaijx

(k)j − ∑n

j=i+1,j 6=i aijx(k)j + bi ⇒

Dx(k+1) = −(L + U)x(k) + b

Gauss-Seidel:∑i−1

j=1,j 6=iaijx

(k+1)j + aiix

(k+1)i = −∑n

j=i+1,j 6=i aijx(k)j + bi ⇒

(L + D)x(k+1) = −Ux(k) + b


Metodi iterativiJacobi e Gauss-Seidel sono casi particolari di “splitting” di A: si cercano duematrici M ,N , con M non singolare tali che

A = M − N

da cuiAx = b ⇔ (M − N)x = b ⇔ Mx = Nx + b ⇔ x = M−1Nx + M−1b

Si ricava il procedimento iterativo

x(k+1) = M−1Nx(k) + M−1b

OSS: Non è detto che per calcolare x(k) si calcoli effettivemante l’inversa diM , in generale si risolve il sistema lineare “facile”

Mx(k+1) = Nx(k) + b

JACOBI: M = D, N = −(L + U);

GAUSS-SEIDEL: M = L + D, N = −U


Metodi iterativi: convergenzaDEF: Un metodo iterativo per sistemi lineari

x(k+1) = Px(k) + q, P : (matrice di iterazione)

e convergente se ∀x(0)

limk→∞

‖e(k)‖ = 0

dove e(k) = x(k) − x∗ e x∗ = Px∗ + q (x∗ soluzione esatta)

e(k) = x(k) − x∗ = Px(k−1) + q − (Px∗ + q) = P (x(k−1) − x∗) = · · ·= P k(x(0) − x∗) = P ke(0)

TEO: e(k) = x(k) − x∗ = P ke(0) converge a 0 ∈ Rn, ∀x(0) se e solo se {P k} e convergente.

COR: {e(k)} converge a 0 se e solo se ρ(P ) < 1.

COR: {e(k)} converge a 0 se vale ‖P‖ < 1 per una qualunque norma naturale ‖ ‖.


Metodi iterativi: convergenza

È importante sapere “quanto velocemente” converge il metodo.Si può tentare di stimare il numero k di iterazioni necessarie per ridurre l’erroreiniziale di un fattore maggiore o uguale a 10−m:

‖e(k)‖ ≤ 10−m‖e(0)‖ ⇔ ‖e(k)‖‖e(0)‖

≤ 10−m

‖e(k)‖ = ‖P ke(0)‖ ≤ ‖P k‖ ‖e(0)‖ per ogni norma naturale indotta dalla normavettoriale.Si può scegliere k tale che ‖P k‖ ‖e(0)‖ ≤ 10−m‖e(0)‖ ovvero ‖P k‖ ≤ 10−m.Tenendo conto che esiste una norma naturale tale che ∀ε > 0

ρ(P k) ≤ ‖P k‖ ≤ ρ(P k) + ε,si cerca k tale che ρ(P k) = ρ(P )k ≤ 10−m ovvero k log10(ρ(P )) ≤ −m da cui

k ≥ m

− log10(ρ(P ))

Si definisce R = − log10(ρ(P )) velocità di convergenza: maggiore è R (più piccolo èil raggio spettrale) minori sono le iterazioni necessarie


Convergenza: condizioni necessarie e sufficienti

Applichiamo i teoremi di convergenza ai metodi di Jacobi e Gauss-Seidel; occorreanalizzare gli autovalori di P = M−1N :

0 = det(λI − M−1N) = det(λM−1M − M−1N) = det(M−1) det(λM − N)

Jacobi:

det

λa11 a12 · · · a1n

a21 λa22

. . . a2n

......

. . . an−1,n

an1 an2 · · · λann

= 0

Gauss-Seidel:

det

λa11 a12 · · · a1n

λa21 λa22

. . . a2n

......

. . . an−1,n

λan1 λan2 · · · λann

= 0

⇒ Non è semplice l’analisi dello spettro di queste matrici.


Convergenza: condizioni sufficienti

Jacobi:

‖M−1N‖∞ = ‖

0 −a12

a11

· · · −a1n

a11

−a21

a22

0. . . −a2n

a22

.... . .

. . ....

− an1

ann− an2

ann· · · 0

‖∞

= maxi

n∑

j=1,j 6=i

∣

∣

∣

∣

aij

aii

∣

∣

∣

∣

< 1,

n∑

j=1,j 6=i

|aij | < |aii|, ∀ i

ovvero A a diagonale dominante in senso forte. Si dimostra anche per latrasposta.TEO: Il metodo di Jacobi converge se A e a diagonale dominante in senso forte oppure se AT

e a diagonale dominante in senso forte

Gauss-Seidel: si dimostra analogamente il seguenteTEO: Il metodo di Gauss-Seidel converge se A e a diagonale dominante in senso forte.

TEO: Il metodo di Gauss-Seidel converge se A e simmetrica e definita positiva.


Convergenza: condizioni necessarie

1. tr(P ) =∑n

i=1 λi: se |tr(P )| > n ⇒ ∃λi t.c. |λi| > 1

2. det(P ) = Πni=1λi: se |det(P )| > 1 ⇒ ∃λi t.c. |λi| > 1

⇒ per avere la convergenza è necessario( ma non sufficiente) che

|tr(P )| < n, |det(P )| < 1


Convergenza: casi particolari

TEO (Stein-Rosemberg)Se aii 6= 0 e siano maggiori o uguali a zero gli elementi di J , matrice di iterazione del metodo

di Jacobi. Indicando con GS la matrice di iterazione del metodo di Gauss-Seidel, vale una e

una sola delle seguenti:

1. ρ(GS) = ρ(J) = 1

2. ρ(GS) = ρ(J) = 0

3. ρ(GS) < ρ(J) < 1

4. ρ(GS) > ρ(J) > 1

Molte classi di matrici note (ad esempio quelle derivanti dalla discretizzazionedi PDE di tipo ellittico) soddisfano le ipotesi del teorema sopra. In tali situazioniconviene utilizzare il metodo di Gauss-Seidel che converge più velocemente.

Nel caso di matrici tridiagonali è possibile stabilire esattamente quanto ilmetodo di Gauss-Seidel converga (in questo caso) più velocemente delmetodo di Jacobi, infatti vale

ρ(GS) = ρ(J)2

e quindi Gauss-Seidel ha velocità doppia.SISTEMI LINEARI QUADRATI:METODI ITERATIVI – p.28/54

Varianti del metodo GS: metodo SORPer aumentare la velocità di convergenza nei metodi visti si può introdurre unparametro ω variabile, da scegliere in maniera opportuna.Dal metodo di Gauss-Seidel si considera:

x(k+1)i =

ω

aii

bi −i−1∑

j=1,j 6=i

aijx(k+1)j −

n∑

j=i+1,j 6=i

aijx(k)j

+(1 − ω)x(k)i aii 6= 0

x(k+1)i è ottenuto come combinazione lineare di x

(k)i e della i-esima

componente del vettore ottenuta con il metodo di Gauss-Seidel.

In termini matriciali:

aiix(k+1)i + ω

i−1∑

j=1,j 6=i

aijx(k+1)j = −ω

n∑

j=i+1,j 6=i

aijx(k)j + (1 − ω)aiix

(k)i + ωbi

(D + ωL)x(k+1) = [−ωU + (1 − ω)D]x(k) + ωb

Mωx(k+1) = Nωx(k) + ωbMω := D + ωL, Nω := −ωU + (1 − ω)D


Metodo SOR

Mωx(k+1) = Nωx(k) + ωbMetodo di rilassamento

per ω = 1 il metodo coincide con il procedimento di Gauss-Seidel e per ω 6= 1

il numero di operazioni richieste per effettuare un iterazione è dello stessoordine di quelle richieste dal metodo di Gauss-Seidel, in particolare per

ω < 1 si parla di sottorilassamento

ω > 1 si parla di sovrarilassamento

La convergenza e la velocità di convergenza del metodo sono legate a

ρ(M−1ω Nω)

che per A data dipende da ω. Per valori di ω ∈ R tali che 0 < ω < 2 èassicurata la convergenza (Teo. di Kahan)

Per la velocità, si può cercare il parametro di rilassamento ottimo, cioè:

ω t.c. ρ(M−1ω Nω) ≤ ρ(M−1

ω Nω), ∀ω

È un complicato problema di minimo, ma per alcune classi di importantiproblemi (dati da discretizzazione di PDE) il parametro di rilassamento è noto.


Stabilità e costo computazionale

I metodi iterativi sono stabili: la convergenza (per ogni vettore iniziale)“corregge” gli errori di arrotondamento

Ad ogni iterazione sono richieste n2 operazioni: affinchè sia vantaggiosorispetto ai metodi diretti è necessario che sia limitato il numero di iterazioninecessarie

Il grande vantaggio risiede nell’occupazione di memoria: per matrici sparsepossono essere memorizzati solo gli elementi diversi da zero.

Schema

1. Dati M−1N, M−1b, Nmax, toll, x(0)

2. Per k = 1, 2, . . . ,Nmax

(a) x(1) := M−1Nx(0) + M−1b

(b) se ‖x(1) − x(0)‖ < toll ‖x(1)‖⇒ STOP criterio di arresto soddisfatto,

la soluzione e x(1)

altrimenti x(0) := x(1), prosegui

3. STOP: Criterio non soddisfatto in Nmax itezioni


ESERCIZI


Esecizio 1.Studiare la convergenza del metodo di Jacobi e del metodo Gauss-Seidel per la risoluzione del sistema

lineare Ax = b con

A =

2 1 1

1 2 1

1 1 2

Sia A = D + L + U :

2 1 1

1 2 1

1 1 2

=

2 0 0

0 2 0

0 0 2

+

0 0 0

1 0 0

1 1 0

+

0 1 1

0 0 1

0 0 0

1) Per Jacobi consideriamo il seguente splitting di A = M − N con M = D e N = −(L + U) e

M−1N =

0 − 12

− 12

− 12

0 − 12

− 12

− 12

0

= −1

2

0 1 1

1 0 1

1 1 0


Esecizio 1 (segue ...)

Si cercano gli autovalori di M−1N :

det

−λ 1 1

1 −λ 1

1 1 −λ

= −λ3 +3λ+2 = λ(1−λ2)+2(λ+1) = (λ+1)(λ−λ2 +2)

dove (λ + 1)(λ − λ2 + 2) = 0 per λ = −1, 2. Quindi gli autovalori di M−1N sono − 12{−1, 2} e

il raggio spettrale

ρ(M−1N) =1

2max{1, 2} = 1 ⇒ non converge

ad esempio

x(0) =

1

1

1

, b =

0

0

0

⇒ x(1) = M−1Nx(0) = −1

2

0 1 1

1 0 1

1 1 0

1

1

1

= −

1

1

1

al passo k−esimo si ha x(k) = (−1)kx(0), quindi non converge.


Esecizio 1 (segue ... )

2) Per Gauss-Seidel consideriamo il seguente splitting di A = M − N con M = D + L e N = −U

per cui

M−1N =

2 0 0

1 2 0

1 1 2

−1

0 −1 −1

0 0 −1

0 0 0

=

12

0 0

− 14

12

0

− 18

− 14

12

0 −1 −1

0 0 −1

0 0 0

=

0 − 12

− 12

0 14

− 14

0 18

38

ρ(M−1N) = det

−λ − 12

− 12

0 14− λ − 1

4

0 18

38− λ

= −λ det

14− λ − 1

418

38− λ

= −λ

[(

1

4− λ

) (

3

8− λ

)

+1

32

]

= −λ

(

λ2 − 5

8λ +

1

8

)


Esecizio 1 (segue ... )

−λ(

λ2 − 58λ + 1

8

)

= 0 per λ1 = 0, λ2,3 = 116

(

5 ± i√

7)

.

Gli autovalori di M−1N sono{

0, 116

(

5 ± i√

7)}

e il suo raggio spettrale

ρ(M−1N) = max

{

0,

∣

∣

∣

∣

1

16

(

5 ± i√

7)

∣

∣

∣

∣

}

=1

16

√25 + 7 =

4

16

√2 =

√2

4< 1

⇒ il metodo di Gauss-Seidel converge.

OSS: Per il metodo Gauss-Seidel la convergenza e garantita se A e simmetrica e definita positiva.

Poiche A e simmetrica verifico che sia definita positiva controllando il segno di tutti i suoi minori

principali:

det(A1) = a11 = 2 > 0, det(A2) = det

2 1

1 2

= 3 > 0

det(A3) = det(A) = 8 + 1 + 1 − 2 − 2 − 2 = 4 > 0


Esecizio 2.Studiare la convergenza del metodo di Jacobi e del metodo Gauss-Seidel per la risoluzione del sistema

lineare Ax = b con

A =

1 2 −2

1 1 1

2 2 1

Sia A = D + L + U :

1 2 −2

1 1 1

2 2 1

=

1 0 0

0 1 0

0 0 1

+

0 0 0

1 0 0

2 2 0

+

0 2 −2

0 0 1

0 0 0

1) Per Jacobi consideriamo il seguente splitting di A = M − N con M = D e N = −(L + U) e

J := M−1N = I ·

0 −2 2

−1 0 −1

−2 −2 0



si cerca lo spettro di J , calcolando det(J − λI):

det

−λ −2 2

−1 −λ −1

−2 −2 −λ

= −λ3 − 4 + 4 − 4λ + 2λ + 2λ = −λ3

−λ3 = 0, per λi = 0, i = 1, 2, 3 ⇒ ρ(J) = 0 < 1 ⇒ il metodo di Jacobi converge

2) Per Gauss-Seidel consideriamo il seguente splitting di A = M − N con M = D + L e N = −U e

GS := M−1N = (D + L)−1(−U) =

1 0 0

1 1 0

2 2 1

−1

0 −2 2

0 0 −1

0 0 0

=

1 0 0

−1 1 0

0 −2 1

0 −2 2

0 0 −1

0 0 0

=

0 −2 2

0 2 −3

0 0 2



calcoliamo lo spettro di GS:

det(GS − λI) =

−λ −2 2

0 2 − λ −3

0 0 2 − λ

= −λ(λ − 2)2

det(GS − λI) = 0 per λ1 = 0, λ2 = λ3 = 2 ⇒ ρ(GS) = 2 > 1

⇒ il metodo di Gauss-Seidel non converge

Verifichiamo se le condizioni necessarie per la convergenza erano soddisfatte:

JACOBI:

|det(J)| = −4 + 4 = 0 < 1 ⇒ condizione soddisfatta

|tr(J)| = 0 < 3 ⇒ condizione soddisfatta

⇒ il metodo puo convergere

GAUSS-SEIDEL:

|det(GS)| = 0 < 1 ⇒ condizione soddisfatta

|tr(GS)| = 4 > 3 ⇒ condizione non soddisfatta

⇒ il metodo non puo convergere


Esecizio 3.Assegnata la matrice A ed il vettore b:

A =

5 2 0

1 −6 1

1 0 4

, b =

12

−2

10

1) si vuol approssimare la soluzione del sistema Ax = b calcolando il termine x(3) della successione

ottenuta con il metodo iterativo di Jacobi a partire dal vettore iniziale x(0) = [1 0 0]T .

1) Il sistema lineare da risolvere e:

5x1 + 2x2 = 12

x1 − 6x2 + x3 = −2

x1 + 4x3 = 10

⇒

x(k+1)1 = (12 − 2x

(k)2 )/5

x(k+1)2 = (−2 − x

(k)1 − x

(k)3 )/(−6)

x(k+1)3 = (10 − x

(k)1 )/4

Partendo dal vettore x(0) = [1 0 0]T si ottengono immediatamente i termini successivi:

x(1) =

[

12

5

1

2

9

4

]T

, x(2) =

[

11

5

133

120

19

10

]T

, x(3) =

[

587

300

61

60

39

20

]T



2) Tenendo conto che la soluzione del sistema e data x = [2 1 2]T , e possibile dare una giustificazione

della buona approssimazione ottenuta con x(3)?

2) Per Jacobi si considera il seguente splitting di A = M − N con M = D e N = −(L + U) e

J := M−1N =

5 0 0

0 −6 0

0 0 4

−1

0 −2 0

−1 0 −1

−1 0 0

=

0 −2/5 0

1/6 0 1/6

−1/4 0 0

Si vede immediatamente che ‖M−1N‖∞ = 25

< 1 (oppure ‖M−1N‖1 = 512

< 1) e quindi il

metodo di Jacobi applicato a questa matrice ha una buona velocita di convergenza. Inoltre, senza

entrare nell’analisi dello spetto della matrice di iterazione, si puo dire che il metodo di Jacobi e

convergente perche la matrice A e a diagonale dominante in senso forte.


Esecizio 4.Assegnata la matrice A e il vettore b:

A =

2 0 1

1 2 0

1 0 3

, b =

0

0

1

;

1) si chiede di approssimare la soluzione del sistema lineare Ax = b calcolando il termine x(3) della

successione ottenuta con il metodo iterativo di Gauss-Seidel a partire dal vettore iniziale

x(0) = [1, 0, 1]T ;

2) determinare inoltre se il metodo converge e in caso affermativo determinarne la velocita di

convergenza.

3) In generale, sotto quali condizioni sulla matrice A e assicurata la convergenza del metodo di

Gauss-Seidel?


Esecizio 4 (segue ...)1) Il sistema lineare da risolvere e il seguente:

2x1 + x3 = 0

x1 + 2x2 = 0

x1 + 3x3 = 1

con soluzione 110

[−2, 1, 4]. Dall’ applicazione delle formule di Gauss-Seidel si ottiene:

x(k+1)1 = −x

(k)3 /2

x(k+1)2 = −x

(k+1)1 /2

x(k+1)3 = (1 − x

(k+1)1 )/3

partendo dal vettore iniziale x(0) = [1, 0, 1]T , si ottengono i termini successivi:

x(1) =

[

−1

2,1

4,1

2

]T

, x(2) =

[

−1

4,1

8,

5

12

]T

, x(3) =

[

− 5

24,

5

48,29

72

]T


Esecizio 4 (segue ...)2) Un metodo iterativo per sistemi lineari puo essere espresso

x(0) ∈ Rn,

x(k+1) = M−1Nx(k) + M−1b .

dove A = M − N , det(M) 6= 0. Si ricorda che il metodo converge qualunque sia il vettore iniziale se e

soltanto se ρ(M−1N) < 1 e che la velocita di convergenza e data da R = − log(ρ(M−1N)). Si

ricorda inoltre che ρ(M−1N) ≤ ||M−1N || per ogni norma matriciale naturale. In particolare nel

metodo di Gauss-Seidel si ha:

M =

2 0 0

1 2 0

1 0 3

, N =

0 0 −1

0 0 0

0 0 0

da cui

M−1N =

12

0 0

− 14

12

0

− 16

0 13

.

0 0 −1

0 0 0

0 0 0

=

0 0 − 12

0 0 14

0 0 16



Poiche

||M−1N ||∞ = max

{

1

2,1

4,1

6

}

=1

2< 1,

il metodo di Gauss-Seidel applicato al sistema lineare dato converge.

Si poteva anche dire che il metodo di Gauss-Seidel e convergente perche la matrice A e a

diagonale dominante in senso forte.

Gli autovalori di M−1N sono {0, 0, 16} per cui

ρ(M−1N) =1

6

e la velocita di convergenza

R = − log10(1/6) = 1.7918

.


Esecizio 5.Assegnata la matrice A e il vettore b:

A =

8 0 1

2 4 1

1 0 2

, b =

9

11

3

;

1) si chiede di approssimare la soluzione del sistema lineare Ax = b calcolando il termine x(3) della

successione ottenuta con i metodo iterativo di Jacobi a partire dal vettore iniziale

x(0) = [0, 1, 0]T ;

2) determinare inoltre se il metodo converge qualunque sia il vettore iniziale e in caso affermativo

determinarne la velocita di convergenza.

1) Il sistema lineare da risolvere e il seguente:

8x1 + x3 = 9

2x1 + 4x2 + x3 = 11

x1 + 2x3 = 3


Esecizio 5 (segue ...)dall’ applicazione delle formule di Jacobi si ottiene:

x(k+1)1 = (9 − x

(k)3 )/8

x(k+1)2 = (11 − 2x

(k)1 − x

(k)3 )/4

x(k+1)3 = (3 − x

(k)1 )/2

partendo dal vettore iniziale x(0) = [0, 1, 0]T , si ottengono i termini successivi:

x(1) =

[

9

8,11

4,3

2

]T

, x(2) =

[

15

16,29

16,15

16

]T

, x(3) =

[

129

128,131

64,33

32

]T

2) Un metodo iterativo per sistemi lineari puo essere espresso

x(0) ∈ Rn,

x(k+1) = M−1Nx(k) + M−1b .

dove A = M − N , det(M) 6= 0. Si ricorda che il metodo converge qualunque sia il vettore

iniziale se e soltanto se ρ(M−1N) < 1 e che la velocita di convergenza e data da

R = − log(ρ(M−1N)).



Si ricorda inoltre che ρ(M−1N) < ||M−1N || per ogni norma matriciale naturale. In particolare

nel metodo di Jacobi si ha M = diag(A), quindi:

M =

8 0 0

0 4 0

0 0 2

, N =

0 0 −1

−2 0 −1

−1 0 0

da cui

M−1N = −1

8

0 0 1

4 0 2

4 0 0

.

Poiche ||M−1N ||∞ = 18

max{1, 6, 4} = 34

< 1, il metodo di Jacobi applicato al sistema

lineare dato converge. Gli autovalori di M−1N sono {0, 14,− 1

4} per cui ρ(M−1N) = 1

4e la

velocita di convergenza R = − log10(1/4) = 0.6021.

Si poteva anche dire che il metodo di Jacobi e convergente perche la matrice A e a diagonale

dominante in senso forte.


Esecizio 6Verificare le ipotesi del terorema di Stein-Rosemberg per la matrice A:

A =

6 0 0

−7 9 0

−4 −1 8

,

e confrontare la convergenza del metodo di Jacobi e del metodo di Gauss-Seidel applicati ad A

analizzando i rispettivi raggi spettrali.

La condizione aii 6= 0 e soddisfatta per i = 1, 2, 3. Calcoliamo la matrice di iterazione di Jacobi.

M =

6 0 0

0 9 0

0 0 8

, N =

0 0 0

7 0 0

4 1 0

J := M−1N =

0 0 0

79

0 0

12

18

0

.

La matrice J e a termini positivi per cui le ipotesi del teoroma sono verificate.


Esecizio 6 (segue ...)Calcoliamo il raggio spettrale di J:

det

−λ 0 0

79

−λ 0

12

18

−λ

= −λ3

per cui ρ(J) = 0, e per il teorema di Stein-Rosemberg anche il raggio spettrale di Gauss-Seidel sara nullo.

Verifichiamo:

M =

6 0 0

−7 9 0

−4 −1 8

, N =

0 0 0

0 0 0

0 0 0

⇒ GS := M−1N =

0 0 0

0 0 0

0 0 0

⇒ ρ(GS) = 0. Quindi covergono entrambi alla massima velocita.



A =

9 −3 −1

−2 9 0

−2 0 9

,




M =

9 0 0

0 9 0

0 0 9

, N =

0 3 1

2 0 0

2 0 0

J := M−1N =

0 13

19

29

0 0

29

0 0

.

La matrice J e a termini positivi per cui le ipotesi del teoroma sono verificate.


Esecizio 7 (segue ...)Calcoliamo il raggio spettrale di J:

det

−λ 13

19

29

−λ 0

29

0 −λ

= −λ3 +2

81λ +

2

27λ = λ

(

−λ2 +8

81

)

per cui gli autovalori sono λ1 = 0, λ2,3 = ± 29

√2 ' ±0.3143 ⇒ ρ(J) ' 0.3143 < 1, e per il

teorema di Stein-Rosemberg anche il raggio spettrale di Gauss-Seidel e minore 1 con

ρ(GS) < ρ(J) < 1. Verifichiamo:

M =

9 0 0

−2 9 0

−2 0 9

, N =

0 3 1

0 0 0

0 0 0

GS := M−1N =

19

0 0

281

19

0

281

0 19

0 3 1

0 0 0

0 0 0

=

0 13

19

0 227

281

0 227

281

=1

81

0 27 9

0 6 2

0 6 2


Esecizio 7 (segue ...)Calcoliamo il raggio spettrale di GS:

det

−λ 27 9

0 6 − λ 2

0 6 2 − λ

= −λ[(6 − λ)(2 − λ) − 12] = −λ2(λ − 8)

per cui gli autovalori sono λ1 = λ2 = 0, λ3 = 8 ⇒ ρ(GS) = 881

' 0.0988.



A =

1 0 0

−8 1 −2

−6 −3 6

,




J =

0 0 0

8 0 2

1 12

0

, det

−λ 0 0

8 −λ 2

1 12

−λ

= −λ3 + λ

per cui gli autovalori sono {0, 1,−1} e ρ(J) = 1 per il terorema di Stein-Rosemberg anche ρ(GS) = 1

quindi i due metodi applicati ad A non convergono


SISTEMI LINEARI QUADRATI: METODI ITERATIVI

Documents

Transcript of SISTEMI LINEARI QUADRATI: METODI ITERATIVI