1

Appunti di Calcolo Numerico Lezioni8-11

Sistemi lineari

L’esigenza di risolvere un sistema lineare nasce da svariate applicazioni pratiche.

Ad esempio, consideriamo il problema di voler costruire una dieta che soddisfi certi

requisiti.

Si supponga di voler preparare una colazione con latte magro, pane, marmellata e

biscotti in modo da ottenere 476 calorie. Per essere equilibrata la dieta deve contenere

10 grammi di proteine, 80 grammi di carboidrati, 10 grammi di grasso.

La seguente tabella riporta il numero di calorie, le proteine (espresse in grammi), i

carboidrati (espressi in grammi) e i grassi (espressi in grammi), fornita da ogni

grammo di alimento.

Latte Magro Pane Marmellata Biscotti

Calorie 0.37 4.08 3.05 4.51

Proteine 0.04 0.09 0.006 0.067

Carboidrati 0.05 0.70 0.58 0.74

Grassi 0.001 0.08 0.002 0.14

Indichiamo con x1 , x2 , x3 , x4 , rispettivamente, il numero di grammi dei quattro

alimenti necessario per soddisfare le richieste della nostra dieta. Il valore cercato è la

soluzione del seguente sistema lineare:

=+++=+++

=+++=+++

1014.0002.008.0001.0

8074.058.070.005.0

10067.0006.009.004.0

47651.405.308.437.0

4321

4321

4321

4321

xxxx

xxxx

xxxx

xxxx

.

2

La soluzione di questo sistema lineare è: x1 = 102.1570, x2 =19.8330, x3 =29.9922

x4 = 58.9373. Ciò significa che la nostra dieta richiede che ci siano 102.1572 grammi

di latte magro, 19.8330 grammi di pane, 29.9922 grammi di marmellata e 58.9373

grammi di biscotti.

Studieremo in seguito che affinchè questo sistema ammetta una ed una sola soluzione

è necessario studiare le proprietà della sua matrice dei coefficienti.

Generalità sui sistemi lineari

Consideriamo un sistema lineare

Ax=b

con )( nmMA ×∈ , mn RbRx ∈∈ ,

Si possono verificare i seguenti casi:

Caso 1) m<n, cioè abbiamo più incognite che relazioni lineari tra di esse. Il sistema

è del tipo

=

<bx

nm

A

e si dice indeterminato. Se indichiamo con k il rango della matrice A e k<n, allora il

sistema ammette kn−∞ soluzioni.

3

2) Caso m>n, cioè abbiamo meno incognite che relazioni lineari tra di esse. Il sistema

è del tipo:

=

>bx

nm

A

e si dice sistema sovradeterminato. In questo caso il sistema non ammette una

soluzione esatta, ma una soluzione approssimata. Esiste la soluzione esatta solo se il

rango(A) è uguale al rango [A | b].

3) Caso m=n, cioè il numero di incognite è uguale al numero di relazioni lineari tra di

esse. Il sistema è della forma

=

=bx

nm

A

e si dice sistema normale e sotto opportune ipotesi può ammettere una ed una sola

soluzione.

D’ora in avanti parleremo solo dei sistemi normali, con n=m.

Teorema: Condizione necessaria e sufficiente affinché il sistema lineare Ax=b ,

)( nnMA ×∈ , nRbx ∈, ammetta una ed una sola soluzione, comunque si scelga b, è

che la matrice A sia a rango massimo (cioè che la matrice A sia invertibile); si ha

perciò:

bAx 1−=

4

Un metodo ben noto per la soluzione di sistemi lineari è basato sulla regola di

Cramer, che calcola la componente j-esima della soluzione facendo uso del calcolo

del determinante (realizzato mediante la formula di Laplace). Questo metodo

comporta una complessità computazionale di circa (n+1)! operazioni. Questo

significa che per risolver un sistema con n=10 occorrerebbe fare 39916800

operazioni. Per n=20 (n+1)!=1021. Quindi se sono necessari 10-9 secondi per fare un

prodotto, per risolvere un sistema di n=20 equazioni in n=20 incognite sono necessari

1012 secondi, che equivalgono a circa 1/3 di 105 anni.

Nella pratica quindi vengono utilizzati altri metodi per la soluzione di un sistema

lineare.

Metodi numerici per la soluzione di un sistema lineare.

I metodi per la soluzione di sistemi lineari vengono usualmente divisi in due

raggruppamenti:

1) Metodi diretti ⇒ Questi metodi, in assenza di errori di arrotondamento,

conducono alla soluzione esatta in un numero finito di passi. Essi sono adatti per la

soluzione di sistemi con matrice dei coefficienti densa e di moderate dimensioni.

2) Metodi iterativi ⇒ Questi metodi generano una successione di soluzioni, che,

sotto opportune ipotesi, convergono alla soluzione del sistema. La matrice dei

coefficienti non viene modificata durante il calcolo e quindi è più agevole sfuttarne la

sparsità. Sono adatti, quindi, per la soluzione di sistemi con matrice dei coefficienti

di grandi dimensioni e sparsa. In assenza di errori di arrotondamento conducono alla

soluzione esatta in un numero infinito di passi.

Poiché l’analisi di entrambe le classi di metodi andrebbe al di là dello scopo di questo

corso, di seguito verrà considerata solo una classe di metodi diretti.

5

Metodi di fattorizzazione.

Questi metodi si basano sul fatto che esistono delle classi di matrici per cui la

soluzione del sistema

Ax=b è immediata.

Consideriamo ad esempio le matrici triangolari.

Sia

=

nnnnn llll

lll

ll

l

L

...

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

0...0

0......0

321

333231

2221

11

una matrice triangolare inferiore ed

=

nn

n

n

n

r

rr

rrr

rrrr

R

...000

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

...00

...0

...

333

22322

1131211

una matrice triangolare superiore

La soluzione di sistemi lineari con queste matrici dei coefficienti si ottiene facilmente

mediante sostituzione.

Soluzione del sistema Lx=b, nel caso in cui la matrice L sia triangolare inferiore

Sostituzione in avanti

=++−=

−=

=

−− nklxlxlxlbx

lxlbx

l

bx

kkkkkkkkk ,...,2/))...((

...

/)(

,11,2211

2212122

11

11

Cioè, in forma algoritmica:

for i=1,2,..,n

xi=bi

for j=1,2,…,i-1

6

xi=xi-l ij⋅xj

end for j

xi=xi/aii

end for i

Soluzione del sistema Rx=b, nel caso in cui la matrice R sia triangolare superiore

Sostituzione all’indietro

−−=++−=

−=

=

++++

−−−−−

1,...,2,1/))...((

...

/)(

,,22,11,

1,1,111

,

nnkrxrxrxrbx

rxrbx

r

bx

kknnkkkkkkkkk

nnnnnnn

nn

nn

cioè, in forma algoritmica,

for i=n,n-1,...,1

xi=bi

for j=i+1,…,n

xi=xi-r ij⋅xj

end for j

xi=xi/aii

end for i

Complessità computazionale:

Calcoliamo la complessità computazionale dell’algoritmo di sostituzione in avanti.

Per calcolare la prima componente della soluzione è necessaria 1 divisione. Per

calcolare la seconda componente della soluzione sono necessarie 1 moltiplicazione ed

1 divisione, quindi 2 operazioni moltiplicative

Per calcolare la componente i-esima della soluzione sono necessarie i-1

moltiplicazioni ed 1 divisione, quindi i operazioni moltiplicative.

Riassumendo:

7

soluzione Operazioni moltiplicative

x1 1

x2 2

… …

xi i

xn n

La complessità computazionale in termini di operazioni moltiplicative è :

∑=

+=n

i

nni

1 2)1( , quindi dell’ordine di

2

2nO .

Lo stesso risultato vale per la complessità computazionale dell’algoritmo di

sostituzione all’indietro.

Una classe di metodi di fattorizzazione consiste nel fattorizzare la matrice A nel

prodotto di due matrici triangolari, una inferiore L ed una superiore R, tali che L ed R

siano a rango massimo, cioè se

A=LR

si può spezzare il problema originale

Ax=b

in due problemi di immediata soluzione

Ax=LRx=b

che si può spezzare nella forma:

8

==

yRx

bLy.

Il primo problema che sorge per ottenere questa fattorizzazione in modo univoco è

che la matrice A ha 2n elementi, L e R hanno 2

)1( +nn elementi ciascuno, cioè nn +2

in totale; si devono, quindi, imporre n condizioni.

I) Se si impone che gli elementi diagonali di R siano tutti uguali ad 1 si parla

di fattorizzazione di Crout;

II) Se si impone che gli elementi diagonali di L siano tutti uguali ad 1 si parla

di fattorizzazione di Doolittle.

Il secondo problema è chiedersi se assegnata A, a rango massimo, esiste ed è unica la

sua fattorizzazione LR.

La risposta a questo problema è data dal seguente teorema:

Esistenza della fattorizzazione.

Teorema 1 Se i minori principali di ordine k, k=1,…,n della matrice )( nnMA ×∈

hanno rango massimo, allora esiste una ed una sola matrice triangolare inferiore L,

con 1 sulla diagonale, una ed una sola matrice triangolare U con 1 sulla diagonale e

una ed una sola matrice diagonale D=diag(di), di≠0, per cui

A=LDU

E ponendo

R=DU

si ha

9

A=LR

Questa fattorizzazione è unica.

Algoritmo di fattorizzazione di Gauss

L’algoritmo di fattorizzazione di Gauss viene chiamato anche metodo di eliminazione

di Gauss in quanto esso consiste nel sottoporre il sistema lineare che dobbiamo

risolvere ad opportune trasformazioni in modo tale da eliminare successivamente le

incognite dalle varie equazioni fino a ridursi ad un sistema triangolare da risolvere

per sostituzione.

Più precisamente, considerata la matrice completa

[ ]

=

nnnnn

n

n

aaa

baaa

baaa

bA

b|...

...| ..............

..| .........

|...

|....

|

21

222221

111211

Al primo passo si cerca di eliminare l’incognita 1x da tutte le equazioni successive

alla prima; questo significa annullare tutti i termini della prima colonna della matrice

sotto 11a ; per far questo si sottrae un multiplo opportuno della prima equazione dalle

restanti equazioni, in modo tale da annullare in esse i coefficienti di 1x . Lavorando

sugli elementi della matrice, si avrà, quindi, per gli elementi della prima colonna:

niamaa iii ,...,211111 =−=

Da qui si ricava che, se si vuole che questi si annullino, si dovrà scegliere

10

11

11 a

am i

i = dove 011 ≠a per ipotesi

Gli elementi della matrice verranno quindi modificati così:

njbmbb

niamaa

iii

jiijij

,..,2

,...,2

11)1(

11)1(

=−=

=−=

Quindi, dopo il primo passo, la matrice completa diventa:

[ ]

=

−−−−−−−−−−

)1()1()1(2

)2(32

)1(2

)1(2

)1(22

11____

1211

)1(

b|...|0

...| ..............|0

..| ... ...|0

|...|0

| ....

|

nnnn

n

n

aa

a

baa

baaa

bA

Queste operazioni sulla matrice si possono esprimere mediante premoltiplicazione di

[ ]bA | per una matrice di tipo S, (vista nella lezione precedente), ma con gli

elementi della prima colonna tutti diversi da zero, così definita:

−

−−

=

1

...0...

1

01

1

1

31

21

1

nm

m

m

L [ ] [ ]bALbA || 1

)1( =⇒

11

Sulla matrice [ ]bALbA |]|[ 1)1( = operiamo ora come secondo passo ulteriori

trasformazioni atte ad annullare tutti i termini sottodiagonali della seconda colonna.

In termini matriciali questo si può effettuare premoltiplicando [ ]bAL |1 per la

matrice

−

−=

10

.........

10

010

1

2

322

nm

mL dove

22

22 a

am i

i = i=3,4,…,n

Si avrà allora:

[ ]

=

)2()2()2(3

)2(3

(2)3

)2(33

)1(2

)1(2

)1(22

111211

12

b|..00

...| .............

| ...00

|... ...0

| .... ....

|

nnnn

n

n

n

aa

baa

baa

baaa

bALL

In generale, quindi, al k-esimo passo, con 11 −≤≤ nk , verrà eliminata l’incognita

kx dalle equazioni che vanno dalla k+1-esima alla n-esima, mediante le seguenti

trasformazioni della matrice completa:

12

−=+=

−=

+==

−=

kikii

kjikijij

kk

ikik

bmbb

nkji

amaa

nkia

am

nk

,...,1,

,...,1

1,...,1

In termini matriciali ciò equivale a premoltiplicare )1(]|[ −kbA per la matrice

−

−=

+

1

......

10

1

...

01

1

,

,1

kn

kk

k

m

m

L.

Per k=n-1 si avrà:

[ ] [ ]yRbALLLL nn ||.... 1221 =−− (1)

Se poniamo

12211 .... LLLLL nn −−

− =

abbiamo che 1−L è una matrice triangolare inferiore avente tutti gli elementi

diagonali uguali ad 1.

13

Quindi la (1) può essere riscritta come:

[ ] [ ] [ ]yRbLALbAL ||| 111 == −−−

dove

bLyybL

LRARAL

=⇒==⇒=

−

−

1

1

Da cui si vede che alla fine del metodo di eliminazione di Gauss non solo abbiamo

ottenuto la fattorizzazione LR della matrice ma abbiamo anche risolto il sistema

bLybLy 1−=⇒= .

Sarà quindi sufficiente applicare l’algoritmo di sostituzione all’indietro per ottenere

la soluzione del sistema.

E’ da notare che applicando l’algoritmo così come lo abbiamo descritto si ottiene al

posto di A la matrice R, ma non viene memorizzata L.

Infatti, se

( ) 11

12

11

112211221

1 ............ −−

−−−−−−−

− ==⇒= nnnnn LLLLLLLLLLLLL

Ora se

−

−−

=

1

...0...

1

01

1

1

31

21

1

nm

m

m

L

14

si ha

=−

1

...0...

1

01

1

1

31

211

1

nm

m

m

L

e quindi risulta che L sarà data da

=

1...

1

01

1

321

3231

21

nnn mmm

mm

m

L

OMM

Si è soliti utilizzare le locazioni di memoria sottodiagonali di A per memorizzare gli

elementi di L diversi da zero, cioè gli 1,...,1 ,..1 , −=+= nknkimik

Nota bene: L’algoritmo di Gauss così come descritto è realizzabile se e solo se

1,...,1,0 −=≠ nkakk . Questo è garantito dalle ipotesi del Teorema 1. Queste ipotesi

sono molto restrittive e difficili da verificare.

In generale una matrice di cui sappiamo solo che è a rango massimo non è

fattorizzabile nel prodotto LR.

Esempio:

=

13

20A è a rango massimo ma non è fattorizzabile perché 011 =a .

Tuttavia si può dimostrare che esistono classi di matrici per le quali la fattorizzazione

LR esiste: esse sono le matrici a diagonale dominante e le matrici definite positive.

15

Il seguente teorema alleggerisce le ipotesi del Teorema 1 e ci fornisce un nuovo

metodo per fattorizzare matrici quadrate qualsiasi purchè a rango massimo.

Teorema 2: Se )( nnMA ×∈ è a rango massimo, esiste sempre una matrice di

permutazione P per cui:

PA=LR

dove L è triangolare inferiore con 1 sulla diagonale e R è triangolare superiore con

0≠iir , i =1,..,n.

Dal punto di vista algoritmo questo porta all’algoritmo di Gauss con pivotaggio.

Al passo k, prima di calcolare il moltiplicatore mik, se 0=kka , si va a cercare nella

colonna k-esima, a partire dalla riga k-esima, la posizione di riga l in cui si trova il

primo elemento diverso da zero.

Se la riga l è diversa dalla riga k-esima, si effettua lo scambio tra la riga k-esima e la

riga l-esima della matrice e del termine noto, e poi si procede secondo lo schema

dell’algoritmo classico.

Questa scelta garantisce l’esistenza della fattorizzazione LR, per la quale è

sufficiente che 1,...,1,0 −=≠ nkakk

Esiste un’altra variante dell’algoritmo di Gauss, con pivotaggio a perno massimo, che

al passo k, prima di calcolare il moltiplicatore mik, va a cercare nella colonna k-esima,

a partire dalla riga k-esima, la posizione di riga l in cui si trova l’elemento di modulo

massimo.

Vedremo che per la stabilità della fattorizazione è meglio utilizzare la variante del

pivotaggio a perno massimo, in modo tale che gli elementi di L risultano tutti minori

od uguali ad 1.

16

Esempio

Risolviamo con il metodo di eliminazione di Gauss il seguente sistema

Ax=b

=+++=+++

=+++=+++

3520104

201063

10432

4

4321

4321

4321

4321

xxxx

xxxx

xxxx

xxxx

Consideriamo

[ ]

=

35 |201041

20|10631

10| 4321

4|1111

| bA

Al passo k=1, si calcolano i moltiplicatori

111

2121 ==

a

am , 1

11

3131 ==

a

am ,

111

4141 ==

a

am

e si sottrae alla seconda riga la prima riga moltiplicata per 21m ; alla terza riga la

prima riga moltiplicata per 31m , alla quarta riga la prima riga moltiplicata per 41m , la

stessa operazione viene effettuata sul termine noto. Si ottiene:

[ ]

=

13 | 19930

16|9520

6| 3210

4|1111

| )1(bA

17

Al passo k=2, si calcolano i moltiplicatori

222

3232 ==

a

am , 3

22

4242 ==

a

am

e si sottrae alla terza riga la seconda riga moltiplicata per 32m ; alla quarta riga la

seconda riga moltiplicata per 42m ; la stessa operazione viene effettuata sul termine

noto. Si ottiene:

[ ]

=

31 | 10300

4|3100

6| 3210

4|1111

| )2(bA

Al passo k=3, si calcola il moltiplicatore

333

4343 ==

a

am

e si sottrae alla quarta riga la terza riga moltiplicata per 43m ; la stessa operazione

viene effettuata sul termine noto. Si ottiene:

[ ]

=

1 | 1000

4|3100

6| 3210

4|1111

| )3(bA

A questo punto la matrice A è stata ridotta in forma triangolare superiore R ed il

termine noto b è diventato la soluzione y del sistema triangolare inferiore Ly=b.

Bisogna adesso risolvere il sistema triangolare superiore

18

Rx=y

Risolviamolo mediante la sostituzione all’indietro:

Ricaviamo dall’ultima equazione 4x

4x =1

Sostituiamo 4x nella terza equazione e ricaviamo 3x

11/)34( 43 =−= xx

Sostituiamo 3x , 4x nella seconda equazione e ricaviamo 2x

11/)3 x26( 432 =−−= xx

Sostituiamo 2x , 3x , 4x nella seconda equazione e ricaviamo 1x

11/) x4( 4321 =−−−= xxx

Complessità computazionale dell’algoritmo di fattorizzazione di Gauss

Per analizzare il numero di operazioni necessarie per risolvere un sistema lineare

Ax=b mediante l’algoritmo di eliminazione di Gauss, consideriamo separatamente il

numero di operazioni necessarie per la fattorizzazione, per la sostituzione in avanti

bLybLy 1−=⇒= e per la sostituzione all’indietro.

Consideriamo le operazioni necessarie alla fattorizzazione LR: al primo passo sono

necessarie n-1 divisioni per calcolare i moltiplicatori .,...,2, nimi = Poi sono necessarie

2)1( −n moltiplicazioni ed 2)1( −n addizioni per modificare la matrice creando i numeri

)2(ija i=2,…,n j=2,…n.

Sintetizzando si ha:

Passo Addizioni Moltiplicazioni Divisioni

19

1 2)1( −n 2)1( −n (n-1)

2 2)2( −n 2)2( −n 3 (n-2)

… … …. ….

… … … ….

n-1 1 1 1

Totale 6

)12)(1( −− nnn 6

)12)(1( −− nnn 2

)1( −nn

dove abbiamo utilizzato il seguente ben noto risultato.

Dato un intero p≥1

∑∑==

++=+=p

j

p

j

pppj

ppj

1

2

1 6

)12)(1(

2

)1(

Quindi si ha che per la fattorizzazione LR sono necessarie

Moltiplicazioni e divisioni: 33)1( 32 nnn ≈−

Addizioni e sottrazioni: 36)12)(1( 3nnnn ≈−−

Per la soluzione dei due sistemi triangolare sono necessarie 2

2n≈ operazioni

moltiplicative e 2

2n≈ addizioni sottrazioni ciascuno. Quindi complessivamente

2n≈ operazioni moltiplicative e 2n≈ addizioni e sottrazioni.

Quindi, la complessità computazionale totale per risolvere Ax=b sarà data da

20

Moltiplicazioni e divisioni: 32

3

3

1

3nn

n ≈+≈

Addizioni e sottrazioni: 32

3

3

1

3nn

n ≈+≈

Stabilità numerica dell’algoritmo di eliminazione di Gauss per la fattorizzazione

LR

Senza perdere in generalità, consideriamo il caso della fattorizzazione

A=LR

e studiamo l’effetto che ha sui risultati il fatto che questa fattorizzazione venga

eseguita operando con i numeri finiti; poiché le operazioni aritmetiche che

compaiono nell’algoritmo di eliminazione di Gauss per calcolare la fattorizzazione

LR di un a matrice A vengono eseguite in aritmetica finita, alla fine dell’algoritmo si

ottengono, anziché i fattori esatti L ed R i fattori L ed R; questi si possono pensare

come

L=L+δL R=R+δR, (2)

cioè dati dai fattori esatti più una piccola perturbazione.

Utilizzando la filosofia dell’analisi all’indietro introdotta da Wilkinson, i fattori L ed

R possono essere pensati come fattorizzazione esatta di una matrice perturbata:

A+δA= L⋅ R (3)

21

Si studia, perciò, da cosa dipenda l’entità della perturbazione δA.

Dalle relazioni (2) e (3) si ottiene

A+δA= (L+δL ) ⋅ (R+δR)

cioè

A+δA= L⋅R+L⋅δR + δL⋅R +δL⋅δR

da cui segue per δA la seguente espressione:

δA=δL⋅R+ L⋅δR +δL⋅δR.

Questa relazione mette in evidenza che la perturbazione δA, non solo dipende dalle

piccole perturbazioni δL e δR, ma è tanto più grande quanto più grandi sono gli

elementi dei fattori L ed R.

Questa osservazione porta al fatto che si definisce la stabilità della fattorizzazione LR

in termini degli elementi di L e di R.

Definizione di stabilità di un algoritmo di fattori zzazione LR

Si dice che un algoritmo produce una fattorizzazione LR numericamente stabile in

senso forte di una matrice A, i cui elementi sono tutti minori od uguali ad 1, se si

possono trovare delle costanti positive a e b, indipendenti dall’ordine e dagli elementi

di A tali che

bral ijij ≤≤

Se le costanti a e b dipendono dall’ordine di A si dice che la fattorizzazione LR è

stabile in senso debole.

22

Se utilizziamo questa definizione per vedere se l’algoritmo di fattorizzazione di

Gauss si vede che applicando la strategia del pivoting con perno massimo per

colonne si ottiene che tutti gli elementi di L sono minori od uguali ad 1

1≤=kk

ikik a

am ( in quanto dopo l’applicazione della strategia del pivot per

colonne a perno massimo, sicuramente ikkk aa ≥ per i = k+1,…,n).

Quindi, la costante a della definizione risulta essere a=1.

Questo spiega perché si sceglie il massimo elemento in modulo e non semplicemente

un elemento diverso da zero.

Per quanto riguarda gli elementi di R si vede che essi possono crescere

esponenzialmente con n.

Infatti al primo passo:

jiijij amaa 11)1( +≤ con 11 ≤im

e quindi

ijij aa max2)1( ≤

Al secondo passo:

ijijjijjiijij aaaaamaa max22 2)1()1(2

)1()1(22

)1()2( ≤≤+≤+≤

23

Al generico passo k si avrà:

ijkk

ij aa max2)( ≤

Al passo finale k=n-1 risulta:

ijnn

ijij aar max2 1)1( −− ≤=

Quindi, se 1||max,

=ijji

a , la costante b della definizione risulta essere b=2n-1.

Quindi l’algoritmo di Gauss non è fortemente stabile perché b dipende da n.

Wilkinson ha trovato una matrice in cui questa maggiorazione viene raggiunta:

−−−−

−−

−

=

11111

11111

10111

10011

10001

A

Operando la fattorizzazione LR si trova che l’elemento rnn è proprio uguale a 12 −n .

Esempio: Consideriamo il caso in cui la matrice di Wilkinson sia di ordine n=8.

L’elemento rnn è uguale a 128. Effettuiamo una perturbazione di 0.5 sull’elemento

rnn , cioè una perturbazione dello 0.5%.

24

=

5.00000

00000

00000

00000

00000

Rδ

Se andiamo ad effettuare il prodotto L⋅(R+δR) otterremo una matrice A che in

posizione (8,8) ha il valore di 1.5. Cioè abbiamo ottenuto una perturbazione del 50%.

Quindi in generale il metodo di Gauss con pivotaggio parziale produce una

fattorizzazione stabile in senso debole.

Per fortuna nella pratica le cose non vanno poi così male:

a) le valutazione empirica del metodo di Gauss ha mostrato che in generale

2,1nrij ≤≈ , cioè gli elementi di R non crescono poi così velocemente.

b) Esistono classi di matrici per cui la fattorizzazione di Gauss è numericamente

stabile in senso forte.

i. Matrici a diagonale dominante

ii. Matrici simmetriche e definite positive.

25

Condizionamento di un sistema lineare

In questa sezione si vuole esaminare come perturbazioni sugli elementi della

matrice A e sugli elementi del termine noto b influenzano la soluzione x del

sistema lineare. Queste perturbazioni sono tipicamente dovute sia agli errori di

approssimazione quando la matrice A ed il termine noto b vengono rappresentati

con numeri finiti, sia al fatto che tutte le operazioni della fattorizzazione e

soluzione dei sistemi triangolari vengono effettuate in aritmetica finita.

Teorema:

Sia una qualunque norma naturale; sia )( nnMA ×∈ a rango massimo e sia δA una

matrice di perturbazione e δb il vettore di perturbazione del termine noto.

Allora se x è la soluzione di Ax=b e δx soddisfa la relazione:

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

si ha

+≤

A

A

b

bAK

x

x δδδ)( (4)

dove 1)( −= AAAK è detto indice di condizionamento o numero di condizione della

matrice A.

L’indice di condizionamento della matrice identità è uguale ad 1, K(I)=1.

In generale l’indice di condizionamento di una matrice è maggiore o uguale ad 1,

K(A)≥1.

Se K(A) >>1, allora la matrice A è mal condizionata ed il problema Ax=b è mal

posto. In questo caso errori anche molto piccoli sui termini della matrice o del

26

termine noto, possono trasformarsi in errori molto elevati sulla soluzione finale (vedi

(4)).

Se K(A)≈1 allora la matrice A è ben condizionata ed il problema Ax=b è bene posto.

In questo caso errori piccoli sui termini della matrice o del termine noto non vengono

amplificati sulla soluzione.

Esempi di matrici mal condizionate.

1) Matrice di Vandermonde: dato un vettore ),...,,( 10 nxxxx = , la matrice di

Vandermonde è una matrice di dimensione (n+1)×(n+1), il cui generico

elemento di posto (i,j) è dato da

j

iij xa )(= i=0,..n j=0,..n

cioè:

=

23

233

22

222

21

211

30

200

)()(1

)()(1

)()(1

)()(1

xxx

xxx

xxx

xxx

A

2) Matrice di Hilbert

La matrice di Hilbert è così definita:

njiji

hij ,...,1,1

1 =−+

=

Nel caso n=4

=

7

1

6

1

5

1

4

16

1

5

1

4

1

3

15

1

4

1

3

1

2

14

1

3

1

2

11

H

27

L’indice di condizionamente della matrice di Hilbert di ordine 4 è

K(H)=1.55*104.

Osservazione: L’indice di condizionamento di una matrice dipende intrinsecamente

dal problema , dalla matrice stessa, cioè non ha nulla a che vedere con la

fattorizzazione stabile o no; esso ci dice come le inevitabili perturbazioni che si

hanno sulla matrice o sul termine noto si ripercuotono sulla soluzione del sistema,

prescindendo da come questa soluzione si ottiene. Una matrice A simmetrica e

definita positiva può essere mal condizionata; in questo caso anche un algoritmo

fortemente stabile può dare risultati completamente falsi .

Vediamo ora con un esempio banale come un alto indice di condizionamento può

falsare i risultati completamente a partire da piccole perturbazioni sui dati originali.

Consideriamo il sistema lineare Ax=b dove

=

=

7.0

1

75

107bA

La soluzione esatta di questo sistema lineare è x=

1.0

0

Questo punto di R2 il punto di intersezione delle due rette

7.0751107 2121 =+=+ xxexx

28

Sia

−=

01.0

01.0bδ una perturbazione del termine noto.

Risolviamo lo stesso sistema lineare con termine noto perturbato:

=+=

69.0

01.1' bbb δ

La soluzione diventa

−=

22.0

17.0'x

Cioè le rette

69.07501.1107 2121 =+=+ xxexx si intersecano nel punto

−=

22.0

17.0'x

29

Calcoliamo l’indice di condizionamento della matrice A:

−−

=

= −

75

107

75

107 1AA

K(A)=17x17=289

=∞

∞

b

bδ0.01 (errore dell’1% sul termine noto).

−=−=

12.0

17.0'xxxδ

7.11.0

17.0 ==∞

∞

x

xδ

Cioè un’errore relativo sui dati dell’ordine dell’1% comporta un errore relativo sulla

soluzione maggiore del 100% (precisamente del 170%).

30

Nel caso in cui si abbia dinanzi un problema mal posto, si possono seguire

le seguenti strade:

1) Cambiare la formulazione del problema, per superare l’ostacolo.

2) Usare la precisione multipla nei calcoli

3) Usare tecniche di regolarizzazione, che sostituiscono al problema di

partenza un problema leggermente modificato, ma che risulta ben

condizionato.

31

Fattorizzazione di Cholesky per matrici simmetriche e definite positive.

Una matrice A di ordine n si dice simmetrica e definita positiva, se nTT RxxAxxeAA ∈≠∀>= 00

Questa condizione è difficile da verificare, esiste una condizione sufficiente affinchè

una matrice simmetrica sia definita positiva, che risulta utile nella pratica.

Teorema:

Sia nxnRA∈ simmetrica: se A ha elementi diagonali positivi ed è a diagonale

strettamente dominante, cioè se :

∑>≠=

n

ijj

ijii aa1

Allora A è definita positiva.

Non vale il viceversa. Cioè possono esistere matrici A definite positive, che non

soddisfano le ipotesi del teorema.

Per la fattorizzazione delle matrici simmetriche e definite positive è stato studiato un

algoritmo di fattorizzazione, detto di Cholesky, che deriva dal seguente teorema:

Teorema di Cholesky.

Sia A una matrice di ordine n simmetrica e definita positiva, allora esiste una matrice

triangolare inferiore L con elementi diagonali positivi, ( ),...10( nil ii => tale che

TLLA ⋅=

32

Partendo dal Teorema, ricaviamo l’algoritmo che ci permette di costruire L.

Sia A una matrice simmetrica e definita, per il teorema si ha:

=⋅

=

=

nn

niii

njijjj

nij

nij

nnninjnn

iiijii

jjjj

nnninjnn

niiiijii

njijjjjj

nij

nij

l

ll

lll

llll

lllll

lllll

llll

lll

ll

l

aaaaa

aaaaa

aaaaa

aaaaa

aaaaa

A

...0...0...00

000

00

0

.........

0

0

00

00

0000

.........

22222

1112111

21

21

21

2221

11

21

21

21

2222221

1112111

OMMMMMM

MM

MMOMMMM

LM

MMOMM

LLL

LLL

OMMMMMM

MMM

OMMMM

MLM

MOMM

ML

LLL

OMMMMMM

MM

MMOMMMM

LM

MMOMM

LLL

LLL

Imponendo l’uguaglianza termine a termine tra gli elementi della matrice A e gli

elementi della matrice prodotto TLLA ⋅= .

Cominciamo con gli elementi della prima colonna:

1111111111 allla =⇒⋅= (si noti che il radicando sarà sicuramente positivo essendo la

matrice A simmetrica definita positiva per ipotesi. Si prende la radice positiva,

essendo, secondo il teorema, gli elementi diagonali positivi).

33

11

111111 a

allla iiii =⇒⋅= i = 2,…n

Adesso procediamo con gli elementi della seconda colonna:

2212222

222

22122 lallla −=⇒+=

22211222222112 /)( lllallllla jjjjjj −=⇒+= j=2,…n

Costruiamo adesso la colonna j-esima

njlal

lllalllla

jk

j

kjjjj

jjjjjjjjjjjjjj

,...1

......

21

1

1,22

22

122

22

1

=∑−=

⇒−−−−=⇒+++=

−

=

−

njilllal

lllllllal

lllllla

jj

j

kjkikijij

jjjjjijijiiijij

jjijjijiiij

,...1/)(

/)...(

...

1

1

1,1,221

221

+=∑−=

⇒−−−−=

⇒+++=

−

=

−−

34

Quindi, la matrice L viene costruita colonna per colonna, mediante il seguente

algoritmo, detto

Algoritmo di Fattorizzazione di Cholesky:

for j=1,…n

∑−=−

=

1

1

2j

kjkjjjj lal

for i=j+1,..n

jj

j

kjkikijij lllal /)

1

1

∑−=−

=

end for i

end for j

Nota: Se la matrice A è definita positiva, come da ipotesi, i radicandi presenti

nell’algoritmo risultano positivi. L’algoritmo di Cholesky è anche un test per

verificare se la matrice A simmetrica è definita positiva. Se durante il corso

dell’algoritmo si verifica un radicando negativo, si può concludere che la matrice

simmetrica data in input non è definita positiva.

Complessità computazionale.

Per calcolare l’elemento ijl sono necessarie j operazioni moltiplicative.

Sulla riga i-esima della matrice L ci sono solo i elementi, essendo la matrice L

triangolare inferiore.

Quindi per calcolare gli elementi della riga i-esima sono necessarie:

22

)1( 2

1

iiii

i

j≈+=∑

= operazioni moltiplicative.

35

Bisogna costruire i righe, i=1,..,n, quindi la complessità computazionale totale sarà

data da:

3

1

2

6

1

2n

in

i≈∑

=

cioè la metà della complessità computazionale dell’Algoritmo di fattorizzazione di

Gauss.

Analisi della stabilità.

Per studiare la stabilita di un algoritmo di fattorizzazione di una matrice, bisogna,

come abbiamo visto, verificare se esistono delle costanti indipendenti dall’ordine

della matrice che limitano la crescita degli elementi dei fattori in cui è stata

fattorizzata. Nel caso della fattorizzazione di Cholesky, basta fare questo studio per

gli elementi della matrice triangolare inferiore.

Dalla relazione:

221,

221,

22

21 ...... iijiijjiiiii lllllla ++++++= +−

si ricava che

221,

21,

22

21

2 ...... iijijiiiiiij lllllal −−−−−−−= +−

Da cui segue che

iiij al ≤2

Risulta quindi

36

iiij al ≤

Ed indicando con lm l’elemento di motdulo massimo tra gli elementi di L ed con am

l’elemento di modulo massimo degli elementi di A, si può concludere che

mm al ≤

e quindi l’algoritmo di Cholesky è stabile in senso forte perché è possibile

maggiorare gli elementi di L con una costante che non dipende dall’ordine della

matrice.

37

Fattorizzazione di Gauss nel caso di matrici tridiagonali.

Sia A una matrice di ordine n tridiagonale non singolare

=

−−

nn

nnn

iii

ab

cab

cab

cab

cab

ca

A

11

333

222

11

0

0

OOO

OO

A può essere fattorizzata, mediante l’algoritmo di Fattorizzazione di Gauss, nel

prodotto di due matrici L (bidiagonale inferiore) ed R (bidiagonale superiore):

=

=

−−−

n

nn

ii

n

n

i

c

c

c

c

c

RL

αα

α

αα

α

ββ

β

ββ

11

33

22

11

1

3

2

0

0

1

10

1

01

1

OO

O

OOO

O

Applichiamo l’algoritmo di fattorizzazione di Gauss alla matrice A

38

Calcoliamo il moltiplicatore 1

22 a

b=β e sottraiamo alla seconda riga la prima

moltiplicata per 2β .

Gli altri elementi sottodiagonali all’elemento 1,1 sono già nulli.

=

−−

nn

nnn

iii

ab

cab

cab

cab

c

ca

A

11

333

22

11

)1(

0

00

OOO

OO

α

dove 1222 ca βα −=

Il secondo passo dell’algoritmo si riduce semplicemente nel modificare la terza riga

sottraendole la seconda riga moltiplicata per 2

32 a

b=β .

La matrice diventa:

=

−−

nn

nnn

iii

ab

cab

cab

c

c

ca

A

11

33

22

11

)2(

0

0

00

OOO

OO

αα

dove 2333 ca βα −=

Al passo i-1 esimo si modificherà la riga i-esima sottraendole la riga i-1 esima

moltiplicata per

39

1−

=i

ii a

bβ

La matrice diventa:

=

−−

−

nn

nnn

ii

i

ab

cab

c

c

c

ca

A

11

33

22

11

)1(

0

0

00

OOO

O

OO

α

αα

dove

1−−= iiii ca βα

Al passo n-1 esimo si modificherà la riga n-esima sottraendole la riga n-1 esima

moltiplicata per

1−

=n

nn a

bβ

La matrice diventa:

40

=

−

−

n

nn

ii

n

c

c

c

c

ca

A

αα

α

αα

0

0

0

00

1

33

22

11

)1(

OOO

O

OO

dove 1−−= nnnn ca βα .

Quindi la matrice A è stata ridotta a forma triangolare superiore e la matrice L sarà

data da:

=

−

1

10

1

01

1

1

3

2

n

n

i

L

ββ

β

ββ

OOO

O

Di seguito viene riportato l’algoritmo di fattorizzazione di Gauss semplificato per il

caso di matrici tridiagonali:

1

1

11

−

−

−=

=

=

iiii

i

ii

ca

b

a

βαα

β

α

i= 2,…,n

La complessità computazionale dell’algoritmo è:

41

n-1 divisioni, n-1 moltiplicazioni ed n-1 somme.

Quindi 3(n-1) operazioni. E’ quindi lineare.

Se si vuole risolvere un sistema lineare con matrice A tridiagonale, dopo averla

fattorizzata in A=LR, bisogna risolvere i due sistemi lineari:

Ly=b

Rx=y

Essendo L, bidiagonale inferiore, l’algoritmo di sostituzione in avanti per la

soluzione del sistema Ly=b diventa:

y1=b1

yi=bi-βi yi-1 i=2,…n

Essendo R, bidiagonale superiore, l’algoritmo di sostituzione all’indietro per la

soluzione del sistema Rx=y diventa:

iiiii

n

nn

xcyx

yx

αα

/)( 1+−=

= i=n-1,n-2,…1

Ciascuno di questi algoritmi ha complessità dell’ordine di O(n).