Post on 30-Mar-2019
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).