Cn2 Sistemi Lineari Sintesi

Matrici e vettori: cosa non possiamo non sapere Fattorizzazioni di matrice e loro applicazioni Metodi iterativi Convergenza dei meto

Metodi Numerici e Statistici per lIngegneria

Sistemi lineari

Silvia Falletta

Dipartimento di Scienze Matematiche, Politecnico di [email protected]

A.A. 2013/2014


Matrici e vettori: cosa non possiamo non sapere

Sia x = (x1, . . . , xn)T Rn. Richiamiamo le seguenti norme per

vettori:

||x ||2 =n

i=1 x2i =

xT x norma euclidea;

||x || = max1in |xi |;||x ||1 =

ni=1 |xi |.

Sia A = (aij)i ,j=1,...,n Rnn. Richiamiamo le seguenti norme permatrici:

||A||2 =(ATA) norma spettrale, dove (B) = maxi |i |,

con i autovalore di B , e` detto raggio spettrale di B ;

||A|| = max1inn

j=1 |aij |;||A||1 = max1jn

ni=1 |aij |.


Proprieta`Per le norme 1, 2 ed valgono le seguenti proprieta`:

||AB || ||A||||B ||||I || = 1, con I matrice identita`||Ax || ||A||||x ||


Una matrice A si dice triangolare superiore se

A =

? ? ? ? ?0 ? ? ? ?0 0 ? ? ?0 0 0 ? ?0 0 0 0 ?

aij = 0, i > j

Una matrice A si dice triangolare inferiore se

A =

? 0 0 0 0? ? 0 0 0? ? ? 0 0? ? ? ? 0? ? ? ? ?

aij = 0, i < j


Una matrice A si dice diagonale se

A =

? 0 0 0 00 ? 0 0 00 0 ? 0 00 0 0 ? 00 0 0 0 ?

aij = 0, i 6= j

Una matrice A si dice tridiagonale se

A =

? ? 0 0 0? ? ? 0 00 ? ? ? 00 0 ? ? ?0 0 0 ? ?

aij = 0, |i j | > 1


Una matrice A si dice a diagonale dominante (per righe) se

|aii | >n

j=1,j 6=i

|aij |,i = 1, . . . , n

Esempio 4 1 12 7 1

3 2 9

Una matrice A si dice a diagonale dominante per colonne se

|ajj | >n

i=1,i 6=j

|aij |,j = 1, . . . , n

Esempio 6 1 12 7 1

3 2 9


Una matrice A si dice simmetrica definita positiva se

1 AT = A

2 xTAx > 0 x 6= oEsempioLa matrice A = BTB , con B Rnn e non singolare, e` simmetricadefinita positiva.

Una matrice simmetrica A e` definita positiva se, e solo se, tutti isuoi autovalori sono reali e positivi.


Risoluzione numerica di sistemi lineari: condizionamento

Consideriamo il sistema lineare

a11x1 + a12x2 + . . .+ a1nxn = b1a21x1 + a22x2 + . . .+ a2nxn = b2

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

che in forma matriciale diventa

a11 a12 . . . a1na21 a22 . . . a2n...

......

...an1 an2 . . . ann

x1x2...xn

=

b1b2...bn

Ax = b,

e risolviamo il seguente problema: dati A Rnn non singolare(det(A) 6= 0) e b Rn, determinare x Rn tale che Ax = b. Taleproblema ammette una ed una sola soluzione.


Esaminiamo dapprima il condizionamento del problema. A talscopo perturbiamo i dati A e b nei dati A = A + A e b = b + bed indichiamo con x = x + x la soluzione in aritmetica esatta delsistema perturbato Ax = b. Quindi confrontiamo lerrore relativodella soluzione ||x ||/||x || con gli errori relativi dei dati ||A||/||A||e ||b||/||b||. Se ||A|| < 1/||A1|| ( A + A e` non singolare), sidimostra che

||x ||||x ||

K (A)

1 K (A) ||A||||A||

( ||A||||A|| +

||b||||b||

)

doveK (A) = ||A||||A1||.

Supponendo inoltre che ||A|| < 1/(2||A1||) si ha||x ||||x || 2K (A)

( ||A||||A|| +

||b||||b||

).


Pertanto K (A) = ||A||||A1|| viene definito come numero dicondizionamento del sistema Ax = b. Per le norme 1,2 ed dimatrice si ha:

K (A) = ||A||||A1|| ||AA1|| = ||I || = 1.

Pertanto, se K (A) 1 il sistema e` ben condizionato; seK (A) >> 1 il sistema potrebbe essere mal condizionato


Un esempio classico di sistema mal condizionato e` quello associatoalla matrice di Hilbert

Hn =

1 12 . . .1n

12

13 . . .

1n+1

......

......

1n

1n+1 . . .

12n1

, K2(Hn) 10n+1

oppure alla matrice di Vandermonde

Vn =

xn11 . . . x1 1

xn12 . . . x2 1...

......

...xn1n . . . xn 1

dove xi 6= xj per i 6= j .


Comandi MATLAB

norm(A) oppure norm(A,2) calcola la norma 2 del vettore odella matrice A;

norm(A,inf) calcola la norma del vettore o della matriceA;

norm(A,1) calcola la norma 1 del vettore o della matrice A;

cond(A) oppure cond(A,2) calcola il numero dicondizionamento in norma 2 del sistema Ax = b;

cond(A,inf) calcola il numero di condizionamento in norma del sistema Ax = b;cond(A,1) calcola il numero di condizionamento in norma 1del sistema Ax = b;

hilb(n) genera la matrice di Hilbert Hn di ordine n;

vander(x) genera la matrice di Vandermonde associata alvettore x .


Come risolvere numericamente Ax = b?Regola di Cramer

Abbiamo bisogno di una routine per calcolare il determinantedi una matrice.

CATTIVA IDEA: perche` richiede (n + 1)! operazioni dimacchina

Ad esempio, se supponiamo di eseguire un flop (floating pointoperation) in 106 secondi

0 2 4 6 8 10 12 14 16 18 201015

1010

105

100

105

1010

size of the system

year

s


Metodi diretti

I metodi numerici per la risoluzione di un sistema lineare sonoessenzialmente di due tipi:

diretti

iterativi.

Consideriamo dapprima i metodi diretti.


Tecnica di sostituzione allindietro (backward substitution)per sistemi triangolari superiore

Dato il sistema triangolare superiore

a11x1 + a12x2+ . . .+ a1nxn = b1. . . . . . . . .an1n1xn1+ an1nxn = bn1

annxn = bn

, aii 6= 0,i .

ricaviamo direttamente lincognita xn dallultima equazione, xn1dalla penultima equazione,..., x1 dalla prima:

xn =bnann

xn1 =bn1an1nxn

an1n1...

x1 =b1

Pnj=2 a1jxj

a11

{

xn =bnann

xi =bi

Pnj=i+1 aijxj

aii

, i = n 1, . . . , 1

Costo computazionale in termini di operazioni aritmetiche (+, ):n2/2.


Comandi MATLAB

A = [...; ...; ...];

b = [.; .; .];

n = length(b);

x = zeros(n,1);

x(n) = b(n)/A(n,n);

for i = n-1:-1:1

s = 0;

for j = i+1:n

s = s+A(i,j)*x(j);

end

x(i) = (b(i)-s)/A(i,i);

end

s = A(i,i+1:n)*x(i+1:n)


Tecnica di sostituzione in avanti (forward substitution) persistemi triangolari inferiore

Dato il sistema triangolare inferiore

a11x1 = b1a21x1+ a22x2 = b2

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

, aii 6= 0,i

ricaviamo direttamente lincognita x1 dalla prima equazione, x2dalla seconda equazione,..., xn dallultima:

x1 =b1a11

x2 =b2a21x1

a22...

xn =bn

Pn1j=1 anjxj

ann

{

x1 =b1a11

xi =bi

Pi1j=1 aijxj

aii

, i = 2, . . . , n

Costo computazionale in termini di operazioni aritmetiche (+, ):n2/2.


Comandi MATLAB...n=length(b);

x=zeros(n,1);

x(1)=b(1)/A(1,1);

for i=2:n

s=A(i,1:i-1)*x(1:i-1);

x(i)=(b(i)-s)/A(i,i);

end


Metodo delle eliminazioni di Gauss

Assegnato un sistema lineare Ax = b di ordine n, il metodo delleeliminazioni di Gauss trasforma in n 1 passi il sistema Ax = b inuno equivalente (ovvero che ammette la stessa soluzione) Ux = b,con U matrice triangolare superiore.Il metodo di Gauss utilizza le seguenti proprieta` dei sistemi lineari:

la soluzione rimane invariata se si scambiano tra loro dueequazioni del sistema;

la soluzione rimane invariata se si sostituisce ad unequazionedel sistema una combinazione lineare dellequazione stessa conunaltra.


Per la descrizione del metodo delle eliminazioni di Gaussconsideriamo il seguente sistema lineare di ordine n = 4:

Ax = b

a11x1 + a12x2 + a13x3 + a14x4 = b1a21x1 + a22x2 + a23x3 + a24x4 = b2a31x1 + a32x2 + a33x3 + a34x4 = b3a41x1 + a42x2 + a43x3 + a44x4 = b4

STEP k = 1Poniamo a

(1)ij := aij e b

(1)i := bi . Supponiamo a

(1)11 6= 0, altrimenti

scambiamo la prima equazione con lequazione k-esima tale che

a(1)k1 6= 0; quindi eliminiamo lincognita x1 nelle equazioni

i = 2, 3, 4. A tal scopo sostituiamo li -esima equazione conlequazione che si ottiene sommando alli -esima stessa la prima

equazione moltiplicata per mi1 := a(1)i1 /a(1)11 :


mi1 (a(1)11 x1+ a

(1)12 x2 + a

(1)13 x3 + a

(1)14 x4) = mi1b

(1)1

a(1)i1 x1+ a

(1)i2 x2 + a

(1)i3 x3 + a

(1)i4 x4 = b

(1)i

a(2)i2 x2 + a

(2)i3 x3 + a

(2)i4 x4 = b

(2)i

dove

a(2)ij = a

(1)ij + mi1a

(1)1j , b

(2)i = b

(1)i + mi1b

(1)1 , i , j = 2, 3, 4.


STEP k = 2

a(1)11 x1 + a

(1)12 x2 + a

(1)13 x3 + a

(1)14 x4 = b

(1)1

a(2)22 x2 + a

(2)23 x3 + a

(2)24 x4 = b

(2)2

a(2)32 x2 + a

(2)33 x3 + a

(2)34 x4 = b

(2)3

a(2)42 x2 + a

(2)43 x3 + a

(2)44 x4 = b

(2)4

Supponiamo a(2)22 6= 0; quindi eliminiamo lincognita x2 nelle

equazioni i = 3, 4. A tal scopo sostituiamo li -esima equazione conlequazione che si ottiene sommando alli -esima stessa la seconda

equazione moltiplicata per mi2 := a(2)i2 /a(2)22 . Otteniamo

a(1)11 x1 + a

(1)12 x2 + a

(1)13 x3 + a

(1)14 x4 = b

(1)1

a(2)22 x2 + a

(2)23 x3 + a

(2)24 x4 = b

(2)2

a(3)33 x3 + a

(3)34 x4 = b

(3)3

a(3)43 x3 + a

(3)44 x4 = b

(3)4

dove

a(3)ij = a

(2)ij + mi2a

(2)2j , b

(3)i = b

(2)i + mi2b

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


STEP k = 3Supponiamo a

(3)33 6= 0; quindi eliminiamo lincognita x3

nellequazione i = 4. A tal scopo sostituiamo li -esima equazionecon lequazione che si ottiene sommando alli -esima stessa la terza

equazione moltiplicata per mi3 := a(3)i3 /a(3)33 . Otteniamo cos` ilseguente sistema triangolare superiore

a(1)11 x1 + a

(1)12 x2 + a

(1)13 x3 + a

(1)14 x4 = b

(1)1

a(2)22 x2 + a

(2)23 x3 + a

(2)24 x4 = b

(2)2

a(3)33 x3 + a

(3)34 x4 = b

(3)3

a(4)44 x4 = b

(4)4

Ux = b

dove

a(4)ij = a

(3)ij + mi3a

(3)3j , b

(4)i = b

(3)i + mi3b

(3)3 , i , j = 4.

Risolviamo infine il sistema Ux = b con la tecnica di sostituzioneallindietro.


Il seguente schema di calcolo riassume la descrizione del metododelle eliminazioni di Gauss:

i = k+1, . . . , n

k = 1, . . . , n 1

mik = a(k)ik /a(k)kka(k+1)ij = a

(k)ij + mika

(k)kj , j = k + 1, . . . , n

b(k+1)i = b

(k)i + mikb

(k)k

xn =

b(n)n

a(n)nn

xk =b(k)k

Pnj=k+1 a

(k)kj

xj

a(k)kk

, k = n 1, . . . , 1

Costo computazionale in termini di operazioni aritmetiche (+, ):O(n3/3).


Confronto tempi di calcolo: Cramer e GEM

0 2 4 6 8 10 12 14 16 18 201015

1010

105

100

105

1010

size of the system

year

s

Cramers ruleGaussian elimination


EsempioApplichiamo il metodo delle eliminazioni di Gauss al sistema:

Ax = b

2 1 1 20 2 0 11 0 2 10 2 1 1

x1x2x3x4

=

0104

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


STEP k = 1Memorizziamo i moltiplicatori mi1 al posto di ai1, i = 2, 3, 4:

A =

2 1 1 20 2 0 1

1/2 1/2 5/2 20 2 1 1

b =

0104

STEP k = 2Memorizziamo i moltiplicatori mi2 al posto di ai2, i = 3, 4:

A =

2 1 1 20 2 0 1

1/2 1/4 5/2 9/40 1 1 2

b =

01

1/43


STEP k = 3Memorizziamo i moltiplicatori mi3 al posto di ai3, i = 4:

A =

2 1 1 20 2 0 1

1/2 1/4 5/2 9/40 1 2/5 29/10

b =

01

1/429/10

Applichiamo la tecnica di sostituzione allindietro:

x4 = (29/10)/(29/10) = 1x3 = (1/4 9/4)/(5/2) = 1x2 = (1 (1))/2 = 1x1 = (2 + 1 1)/2 = 1


Comandi MATLAB...for k=1:n-1

for i=k+1:n

A(i,k)=-A(i,k)/A(k,k);

for j=k+1:n

A(i,j)=A(i,j)+A(i,k)*A(k,j);

end

b(i)=b(i)+A(i,k)*b(k);

end

end

oppure, piu` semplicemente ed efficientemente,for k=1:n-1

A(k+1:n,k)=-A(k+1:n,k)/A(k,k);

A(k+1:n,k+1:n)=A(k+1:n,k+1:n)+A(k+1:n,k)*A(k,k+1:n);

b(k+1:n)=b(k+1:n)+A(k+1:n,k)*b(k);

end


Pivoting parziale

Nel metodo di Gauss la condizione a(k)kk 6= 0 e` necessaria se si vuole

procedere con le successive eliminazioni. Tale condizione e`garantita se A soddisfa una delle seguenti proprieta`:

det(Ak) 6= 0 k = 1, . . . , n, con Ak matrice di ordine kformata dagli elementi aij , 1 i , j k;A e` a diagonale dominante per righe;

A e` a diagonale dominante per colonne;

A e` simmetrica definita positiva.

Ricordiamo che se a(k)kk = 0 allora occorre individuare lelemento

a(k)ik 6= 0 con k < i n e scambiare la k-esima equazione con la

i -esima. Osserviamo che se A e` non singolare allora esiste

certamente un valore di i per il quale a(k)ik 6= 0.


Per garantire una migliore stabilita` numerica dellalgoritmo diGauss conviene operare uno scambio di equazioni anche quando

|a(k)kk | e` piccolo. Precisamente ad ogni passo k convieneindividuare lindice di riga r per il quale risulta

|a(k)rk | = maxkin

|a(k)ik |

e scambiare la k-esima equazione con lr -esima equazione. Talestrategia e` nota sotto il nome di pivoting parziale. Il pivoting e`superfluo quando:

A e` a diagonale dominante per colonne;

A e` simmetrica a diagonale dominante;

A e` simmetrica definita positiva.


Esempio Consideriamo il sistema di ordine n = 18:

Ax = b, aij = cos((j 1)i ), i = 2i 12n

pi

bi =

nj=1

aij

la cui soluzione esatta e` x = (1, 1, . . . , 1)T . Per esso risulta

K(A) 17||x x ||||x || 10

9

con x calcolata mediante il metodo di eliminazione di Gauss senzapivoting, e

||x xp||||x || 10

15

con xp calcolata mediante il metodo di eliminazione di Gauss conpivoting.


Fattorizzazioni di matrice e loro applicazioni

Il metodo di eliminazione di Gauss con pivoting parziale puo` essereinterpretato come una successione finita di trasformazioni dellamatrice A e del vettore termine noto b. Infatti, al passo k loscambio delle equazioni k ed r del sistema A(k)x = b(k) si puo`realizzare anche nel seguente modo:

PkA(k)x = Pkb

(k)

dovek r

Pk =

1 0 0 0 00 0 0 1 00 0 1 0 00 1 0 0 00 0 0 0 1

k

lr


Inoltre, al passo k la sostituzione delle equazioni i = k + 1, . . . , ncon le equazioni che si ottengono rispettivamente sommando alleequazioni stesse la k-esima equazione moltiplicata per mik , si puo`realizzare nel seguente modo:

MkPkA(k)x = MkPkb

(k)

dovek

Mk =

1 0 0 0 00 1 0 0 00 mk+1k 1 0 00 0 1 00 mnk 0 0 1


Pertanto il metodo delle eliminazioni di Gauss puo` essere cos`ridefinito:

Mn1Pn1 . . .M2P2M1P1 G

Ax =Mn1Pn1 . . .M2P2M1P1 G

b

m

GAx = Gb Ux = bdove

G = Mn1Pn1 . . .M2P2M1P1

eGA = U

prende il nome di fattorizzazione di Gauss.


Riordinando opportunamente i fattori della decomposizioneGA = U, otteniamo

Mn1 . . . M2M1 M

Pn1 . . .P2P1 P

A = U MPA = U

ove le matrici Mk ed Mk sono esattamente dello stesso tipo(triangolare inferiore, con elementi sottodiagonali tutti nulli trannequelli della colonna k dei moltiplicatori e con diagonale unitaria),ma differiscono per un diverso ordinamento dei moltiplicatori mik .Posto

P = Pn1 . . .P2P1

matrice di permutazione,

L = M1


si ha la seguente fattorizzazione:

PA = LU

con L matrice triangolare inferiore con diagonale unitaria ed Umatrice triangolare superiore della decomposizione di GaussGA = U. I fattori della decomposizione PA = LU si ottengonomediante lalgoritmo di Gauss con pivoting parziale e con un costocomputazionale pari a O(n3/3).


Alcune applicazioni della fattorizzazione PA = LU:

risoluzione del sistema lineare Ax = b:

Ax = b PAx = Pb L Uxy

= Pb

{

Ly = Pb yUx = y x

calcolo del determinante di A:

det(A) = (1)sn

i=1

uii ,

ove s e` il numero totale degli scambi di equazioni effettuati;

calcolo dellinversa di A (costo computazionale pari a O(n3)):

PA = LU (PA)1 = (LU)1 A1P1 = U1L1 A1 = U1L1P


risoluzione di p sistemi lineari aventi tutti la stessa matrice deicoefficienti (costo computazionale pari a O(n3/3 + pn2)). Peresempio, per p = 3

Ax1 = b1Ax2 = b2Ax3 = b3

determiniamo PA = LU e quindi calcoliamo xi , i = 1, 2, 3,risolvendo i seguenti due sistemi triangolari:{

Ly = Pbi yUxi = y xi

risoluzione del sistema Apx = b a partire dai dati A e b (costocomputazionale pari a O(n3/3 + pn2)). Per p = 3

A3x = b A (A (Ax)y

z

) = b

Az = bAy = zAx = y


Se A e` simmetrica definita positiva, tenendo conto che il pivoting e`superfluo (P = I ), e che A e` simmetrica, possiamo riscriverePA = LU nella forma

A = LLT

ove L e` una matrice triangolare inferiore con elementi positivi sulladiagonale principale. Tale fattorizzazione e` dettafattorizzazione di Choleski. I fattori della decomposizioneA = LLT si ottengono mediante un algoritmo il cui costocomputazionale e` O(n3/6).Alcune applicazioni della fattorizzazione di Choleski A = LLT :

risoluzione del sistema lineare Ax = b:

Ax = b L LT xy

= b {

Ly = b yLT x = y x

calcolo dellinversa di A (costo computazionale O(2n3/3)):

A = LLT A1 = (LLT )1 A1 = (LT )1L1 A1 = (L1)TL1


Comandi MATLAB

x=A\b calcola x soluzione di Ax = b con il metodo dieliminazione di Gauss con pivoting parziale;

[L,U,P]=lu(A) calcola i fattori L, U, P della fattorizzazionePA = LU di A;

R=chol(A) alcola il fattore R = LT triangolare superiore dellafattorizzazione di Choleski A = RTR = LLT , della matricesimmetrica definita positiva A.


Conclusioni

Di seguito riassumiamo le principali proprieta` dei metodi diretti perla risoluzione del sistema lineare Ax = b:

1 x viene determinata mediante un numero finito di passi (al piu`n 1 se n e` lordine del sistema);

2 x in aritmetica con precisione infinita di calcolo vienedeterminata in maniera esatta; in aritmetica con precisionefinita di calcolo viene determinata con una precisione la cuientita` non dipende dalle richieste dellutente;

3 i metodi diretti modificano la matrice dei coefficienti;

4 i metodi diretti sono efficienti per matrici dense e dipiccole-medie dimensioni.


Metodi iterativi

Consideriamo il sistema lineare

Ax = b

a11x1 + a12x2 + . . .+ a1nxn = b1a21x1 + a22x2 + . . .+ a2nxn = b2

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

e supponiamo che aii 6= 0 i . A partire da un arbitrario vettoreiniziale x(0), con un metodo iterativo si determina una sequenza divettori

{x(k+1)

}k=0,1,...

che, sotto opportune condizioni, convergealla soluzione esatta x del sistema:

x(0) {

x(k)}

k=1,2,...: lim

kx x(k) = o (null vector)


Definiamo due classici metodi iterativi. A tal scopo ricaviamodalli -esima equazione per i = 1, . . . , n la variabile xi :

xi =bi

i1j=1 aijxj

nj=i+1 aijxj

aii, aii 6= 0

e generiamo la seguente successione di vettori:

x(k+1)i =

bi i1

j=1 aijx(k)j

nj=i+1 aijx

(k)j

aii, k = 0, 1, . . .

Tale espressione definisce il cosiddetto metodo di Jacobi.Ponendo invece

x(k+1)i =

bi i1

j=1 aijx(k+1)j

nj=i+1 aijx

(k)j

aii, k = 0, 1, . . .

otteniamo il cosiddetto metodo di Gauss-Seidel.


Esempio{3x1 x2 = 2x1 + 2x2 = 3

={

x1 =2+x23

x2 =3x12

x = (1, 1)T

METODO DI JACOBI:

x

(k+1)1 =

2+x(k)2

3

x(k+1)2 =

3x(k)1

2

x(0) x(1) x(2) x(3) x(4)

0 2/3 7/6 19/18 . . .0 3/2 7/6 11/12 . . .

METODO DI GAUSS-SEIDEL:

x

(k+1)1 =

2+x(k)2

3

x(k+1)2 =

3x(k+1)12

x(0) x(1) x(2) x(3) x(4)

0 2/3 19/18 107/108 . . .0 7/6 35/36 217/216 . . .


I metodi di Jacobi e Gauss-Seidel fanno parte di quella classe piu`generale di procedimenti iterativi che si deducono nel seguentemodo:

Ax = b (D + C )x = b Dx = b Cx= Dx(k+1) = b Cx(k), k = 0, 1, . . .

ove, ad ogni iterazione k, x(k+1) e` la soluzione di un sistemalineare con matrice dei coefficienti D e termine noto b Cx(k).Pertanto, la matrice D deve essere:

1 non singolare, ovvero tale da garantire lesistenza e lunicita`della soluzione ad ogni passo k;

2 di forma semplice (triangolare, diagonale), ovvero tale che lasoluzione del sistema possa ottenersi con un algoritmorelativamente semplice e poco costoso;

3 tale da garantire la convergenza di x(k+1) ad x per k .


Osserviamo che, nel caso del metodo di Jacobi

aiix(k+1)i = bi

i1j=1

aijx(k)j

nj=i+1

aijx(k)j , i = 1, . . . , n,

le matrici D e C sono cos` definite

D =

a11 0 . . . 00 a22 . . . 0...

... . . ....

0 0 . . . ann

, C = A D;


nel caso del metodo di Gauss-Seidel

ij=1

aijx(k+1)j = bi

nj=i+1

aijx(k)j , i = 1, . . . , n

si ha

D =

a11 0 . . . 0a21 a22 . . . 0...

... . . ....

an1 an2 . . . ann

, C = A D


Quando si implementa un procedimento iterativo occorre definiredei criteri darresto. Fissare un numero massimo kmax di iterazionirappresenta senzaltro un criterio darresto.Un altro criterio darresto, che riguarda la bonta`dellapprossimazione x(k), consiste nel fissare una tolleranzarelativa tollr ( eps) o assoluta tolla e nellarrestare il processoalliterazione k se x(k) soddisfa la seguente disuguaglianza

||x(k+1) x(k)|| tollr ||x(k+1)||

oppure||x(k+1) x(k)|| tolla

Osserviamo che tale criterio darresto ha senso quando il metodoconverge rapidamente.


Comandi MATLAB

METODO DI JACOBI:A = [...; ...; ...];

b = [.; .; .];

kmax = . ; toll = .;

n = length(b);

x0 = zeros(n,1);

D = diag(diag(A));

C = A-D;

for k = 1:kmax

x1 = D\(b-C*x0);if norm(x1-x0)


Convergenza dei metodi iterativi

Poniamo e(k) = x x(k) e sottraiamo membro a membro leseguenti due equazioni:

Dx = b CxDx(k) = b Cx(k1)

De(k) = Ce(k1) = e(k) = D1Ce(k1) = e(k) = Be(k1)

doveB = D1C = D1(A D) = I D1A

e` detta matrice di iterazione. Tenendo conto che luguaglianzae(k) = Be(k1) vale per ogni k = 1, 2, . . ., si ha

e(k) = Be(k1) = B2e(k2) = . . . = Bke(0)


da cui deduciamo che

limk

e(k) = o (vettore nullo) x(0) limk

Bk = O (matrice nulla)

Poiche si dimostra che

limk

Bk = O (B) < 1, (B) = maxi|i |, con i autovalore di B

alloralim

ke(k) = o x(0) (I D1A) < 1

Ricordando che per le norme 1,2, sussiste la proprieta`(A) ||A||, possiamo senzaltro affermare che

||I D1A|| < 1 = limk

e(k) = o x(0)


Convergenza dei metodi di Jacobi e Gauss-Seidel

Usando ||I D1A|| < 1 si dimostra che se A e` a diagonaledominante per righe o per colonne, allora i metodi di Jacobi edi Gauss-Seidel convergono.

Si dimostra inoltre che se A e` simmetrica definita positiva,allora il metodo di Gauss-Seidel converge.

Osserviamo infine che la convergenza del metodo di Gauss-Seidelnon implica la convergenza del metodo di Jacobi e viceversa. Seentrambi convergono, in generale il metodo di Gauss-Seidelconverge piu` rapidamente.La rapidita` di convergenza dipende dal raggio spettrale: quanto piu`esso e` piccolo, tanto piu` e` rapida la convergenza.


Conclusioni

Di seguito riassumiamo le proprieta` dei metodi iterativi per larisoluzione del sistema lineare Ax = b:

1 x viene determinata come limite di una successione di vettoriconvergente;

2 x in aritmetica con precisione infinita di calcolo vienedeterminata in maniera approssimata; in aritmetica conprecisione finita di calcolo puo` essere determinata con unaprecisione soddisfacente le richieste dellutente;

3 i metodi iterativi non modificano la matrice dei coefficienti A;

4 i metodi iterativi sono efficienti per matrici sparse e di grandidimensioni.

Matrici e vettori: cosa non possiamo non sapereRisoluzione numerica di sistemi lineari: condizionamentoMetodi direttiMetodo delle eliminazioni di Gauss

Fattorizzazioni di matrice e loro applicazioniMetodi iterativiConvergenza dei metodi iterativi

Cn2 Sistemi Lineari Sintesi

Documents

Transcript of Cn2 Sistemi Lineari Sintesi