La propagazione degli errori nei prodotti tra matrici. · La propagazione degli errori nei prodotti...
Transcript of La propagazione degli errori nei prodotti tra matrici. · La propagazione degli errori nei prodotti...
La propagazione degli errori nei prodotti tra matrici.
Vediamo, in breve, come l’errore fatto ad ogni passo applicando la funzione floating in
prodotti tra matrici incida sul risultato finale. Cominciamo con l’esaminare il prodotto scalare tra
due vettori x, y𝜖 Rn. Se denotiamo con s la variabile destinata a contenere il prodotto scalare
( x, y) = yTx = x1y1 + x2y2 + …xnyn,
questa viene calcolate nel seguente modo. Si calcola s = x1y1, commettendo un primo errore di
arrotondamento al quale ne viene aggiunto un altro, aggiornando s + x2y2 → s. Col successivo
aggiornamento , s + x3y3 →s, si aggiunge un altro errore e così via fino ad avere la variabile finale
s + xnyn →s. In realtà, allora, non viene calcolato il valore vero s, ma un certo numero floating 𝑠 . A
conti fatti si trova la maggiorazione
|𝑠 - s | ≤ nu |y|T |x|+ O(u2),
dove u è l’unità di arrotondamento dell’aritmetica finita. Il termine di ordine superiore O(u2), può
essere eliminato dalla disuguaglianza, supponendo nu < .01, e, modificando in
|𝑠 - s | ≤ 1.01nu |y|T |x|
L’ipotesi nu < .01 è plausibilissima perché se, per esempio, come in ambiente Matlab, u = O(10-16),
il valore .01 è raggiunto se si considerano vettori di lunghezza pari a 10-14 !. Il fatto che a destra ci
siano le somme dei valori assoluti di xiyi, nasce dal fatto che ogni volta che si esegue una
operazione floating, e si calcola fl(a) = a( 1+ ε), si ha fl(a) –a = a ε, da cui |fl(a) –a | ≤ |a| | ε| < |a|
u.
Il fattore nu ( o il fattore 1.01nu), è accettabilissimo, dal momento che si eseguono 2n operazioni
floating e ne risulta un’amplificazione dell’errore di n volte ( non 2n volte come potrebbe
sembrare). Il fattore 2n si ha quando si suppone che le componenti del vettore non siano numeri
floating ma devono essere forzati ad esserlo. In ogni caso che sia n o 2n questo accumulo di errore
è prevedibile ed inevitabile. Nella disuguaglianza, però, interviene un altro fattore che è |y|T|x|,
cioè il prodotto scalare tra i valori assoluti dei due vettori . Questo è un fattore di amplificazione
che merita una certa attenzione. Se passiamo all’errore relativo si trova
𝑠 −𝑠
𝑠 ≤nu
|𝐲|T |𝐱|
|𝐲T𝐱| + O(u2),
oppure
𝑠 −𝑠
𝑠 ≤1.01nu
|y|T |𝐱|
|𝐲T𝐱| ,
dove, il valore s a denominatore a sinistra, viene rimpiazzato da yTx = s, per metterne in luce la
relazione con |y|T|x|. Allora, se accade che |y|T|x| >> |yTx|, il rapporto tra |y|T|x| e |yTx|,
diventa un ampio fattore di amplificazione per l’errore relativo, portando ad un ‘approssimazione
scadente di s.
Esempio ( cancellazione numerica di cifre significative comuni).
Supponiamo che in un algoritmo venga richiesta, ad un certo passo, la differenza tra
numeri molto vicini. Tali numeri siano a e b. La differenza c = a –b, volendo, può essere
interpretata come il prodotto scalare tra il vettore ( 1 1)T e il vettore (a b). Supponiamo che
l’aritmetica finita sia individuata da S = (10,10,?,?), non interessando il più piccolo esponente e il
più grande esponente e che i numeri siano a = 1.234265558967778, b = 1.234265551110447,
possedenti 15 cifre decimali ( o 16 cifre significative), con restanti cifre tutte nulle. Il valore vero è
c = a – b =0.000000007857331 = . 7857331 10-8, messo in floating point standard. Se usiamo
l’aritmetica finita in questione, occorre preliminarmente forzare a e b ad appartenere all’insieme
dei numeri floating. Se l’aritmetica è troncata si ha fl(a) = .1234265558 101 , fl(b) = .1234265551
101, in rappresentazione floating point standard. L’unità di arrotondamento è 10-9 e fl(a) = a( 1 +
εa), fl(b) = b( 1 + εb), | εa| < 10-9, |εb | < 10-9. Il valore effettivamente calcolato è capp =
.0000000007 101 =.7 10-8, dove dopo la cifra 7 le restanti cifre sono tutte nulle. Quindi
rispetto al valore vero l’approssimazione capp cattura solo una cifra significativa !. Ciò è anche
evidenziato dal calcolo dell’errore relativo capp −c
c=-.1091…=O(10-1). In un algoritmo di calcolo il
fatto che un numero molto piccolo sia male approssimato potrebbe avere anche ripercussioni non
negative. Supponiamo , per esempio, che l’algoritmo richieda x = .1 + (a –b). Il valore vero di x è x
= .1 + (a-b) = .100000000785733. Il valore effettivamente calcolato è xapp = .1 + (fla) + fl(b) =
.1000000007. Se confrontiamo i due numeri ci accorgiamo che l’approssimazione cattura ben 10
cifre significative !. In sostanza se un numero piccolo, calcolato con scarsa accuratezza, va a
sommarsi ad un numero grande, calcolata con ottima accuratezza, va a modificare la coda di
questo numero ( le ultime cifre significative) utilizzando la sua testa, cioè le prime cifre
sicuramente corrette. Nell’esempio proposto è bastata la prima cifra 7, che à la stessa del valore
vero, a dare accuratezza alla somma. Ben diversa è la situazione se si deve fare il rapporto di due
numeri piccoli se uno ( o entrambi ) non è accurato.
Supponiamo che il nostro algoritmo richieda y = c/(2 10-8). Il valore vero è y = 7857331/2
=.39286655. Il valore approssimato, invece, è yapp = c(app)/ (2 10-8) = .35. Appena una cifra
significativa corretta !
•
Siano, ora, A e B due matrici e sia n la dimensione comune ( colonne di A, righe di B). Poiché si
tratta di eseguire dei prodotti scalari, il prodotto C = AB è approssimato, in aritmetica finita da C ,
con errore valutato da
C − C ≤ n 𝐮 A B + O(𝐮2)
Anche in questo caso, se vogliamo togliere il pezzo O(u2), con l’ipotesi nu ≤ .01, si può scrivere
C − C ≤ 1.01n 𝐮 A B
Ci aspettiamo un grosso errore relativo se |A | |B | >> |A B | almeno in alcuni coefficienti.. Per
studiare la cosa in modo globale, facciamo riferimento alla norma uno e alla norma infinito. Per
queste norme risulta || |Z| || = ||Z ||, per ogni matrice Z, e vale anche la proprietà di monotonia.
Vale a dire O ≤ X ≤ Y implica ||X | | ≤ ||Y ||.
Segue che
| C − C |
| C |≤ n𝐮
A |B|
| AB |+ O(𝐮2)
Ora la quarta proprietà delle norme ||XY||≤ ||X|| ||Y||, è tale che la disuguaglianza è
abbastanza stretta se X e Y sono matrici non negative, non dovendosi fare somme e differenze
insieme. Allora
| C − C |
| C |≤ n𝐮
A | B |
| AB |+ O(𝐮2)
Allora il fattore di amplificazione dell’errore relativo è il rapporto tra ||A||||B|| e ||AB||. Se
||A||B|| >> ||AB||, ci aspettiamo un grosso errore relativo. Si osservi che si è deciso di scegliere
come errore relativo il rapporto di due norme, cioè misure globali come si è detto in precedenza.
Equazioni lineari.
Sia A una matrice m × n e b un vettore a m componenti. Un classico problema consiste nel
cercare un eventuale vettore x tale che
Ax = b
Più esplicitamente si tratta di un sistema lineare di m equazioni in n incognite in cui la matrice A è la
matrice dei coefficienti e il vettore b è il vettore dei termini noti. 𝑎11𝑥1 + ⋯ 𝑎1𝑗𝑥𝑗 + ⋯ 𝑎1𝑛𝑥𝑛 = 𝑏1
.
.
. 𝑎𝑖1𝑥1 + ⋯ 𝑎𝑖𝑗 𝑥𝑗 + ⋯ 𝑎𝑖𝑛𝑥𝑛 = 𝑏𝑖
.
.
. 𝑎𝑚1𝑥1 + ⋯ 𝑎𝑚𝑗 𝑥𝑗 + ⋯ 𝑎𝑚𝑛 𝑥𝑛 = 𝑏𝑚
Questa è la forma più esplicita di sistema lineare di equazioni. Una diversa utile forma, consiste nello scrivere la matrice A enfatizzando i vettori colonna. Se A la si vede come un susseguirsi di vettori colonna uno dopo l’altro a partire dal primo,
A = (a1 a2 …an),
allora si riconosce subito che il prodotto di A per un certo vettore x a tante componenti quante sono
le colonne di A, consiste nell’eseguire una combinazione lineare delle colonne di A, i cui coefficienti
sono le componenti del vettore. Pertanto, risolvere il sistema Ax = b, significa trovare x a n
componenti tale che
x1a1 + …xnan = b.
Dalle conoscenze sugli spazi lineari di dimensione finita, sappiamo che se i vettori colonna sono
linearmente indipendenti, l’eventuale combinazione lineare che costruisca il vettore b è unica.
Ricordiamo che la lineare indipendenza vuol dire che il vettore nullo si può formare esclusivamente
con coefficienti tutti nulli. O, che è lo stesso, il sistema Ax = 0, ammette solo la soluzione x = 0. Ma
quando è possibile che i vettori colonna sono linearmente indipendenti? Sempre dalle conoscenze
sugli spazi lineari, il numero massimo di vettori in Rn linearmente indipendenti è esattamente n. Un
tipico esempio è quello dei versori e1, …en. Segue che non può essere m < n, ed è possibile che i
vettori colonna di A siano linearmente indipendenti solo se m ≥ n. Purtroppo non possiamo
illustrare il caso m > n in cui si introduce la soluzione nel senso dei minimi quadrati, con svariate
applicazioni in tanti settori dell’ingegneria, della statistica, ecc. Ci limitiamo, dunque, al caso m = n,
supponendo che i vettori colonna di A siano linearmente indipendenti. Come è noto, ogni vettore in
Rn, è generato ( in modo unico) da un sistema di n vettori linearmente indipendenti, da cui il sistema
Ax = b, ammette sempre una sola soluzione. La lineare indipendenza dei vettori colonna di A, si
riconosce dal calcolo del determinante di A. Se questo è nullo i vettori colonna sono linearmente
dipendenti e viceversa. Se det A ≠ 0, il sistema ammette una sola soluzione e viceversa. La matrice
A è detta singolare se detA = 0 e non singolare in caso contrario. Il calcolo di un determinante , di
solito, per piccole dimensioni si effettua in maniera ricorsiva, tenendo presente che il determinante di
uno scalare è lo scalare stesso. Allora se Aij è la matrice ottenuta da A eliminando riga i e colonna j,
si ha
n
1j
ijij
jiAa1A det)(det
Si fissa i e la somma si esegue al variare di J; ma si può fissare j e la somma si esegue al variare di i.
Questa regola non è più praticabile quando la dimensione n comincia a crescere. Infatti se
denotiamo con Pn, il numero di operazioni per il calcolo di un determinante, risulta dalla formula Pn
> Pn-1 ( trascurando le somme). quindi
Pn > nPn-1>n(n-1)Pn-2>…..n(n-1)….3P2,
e, poiché, il calcolo di un determinante di una 2 2 matrice richiede due moltiplicazioni, si ha Pn >
n! Se, per esempio, n = 30, su devono eseguire almeno 30! Operazioni e P30 > 2 × 1032. Se una
moltiplicazione richiedesse 10-9 secondi, occorrerebbero almeno 2×1023 secondi di calcolo !.
Equazioni matriciali.
Nei sistemi lineari di medie dimensioni che illustreremo in queste note, capita molto spesso
di dover trovare la soluzione di molti sistemi lineari con la stessa matrice dei coefficienti. In
sostanza questo vuol dire prendere in esame equazioni matriciali, come, ora, illustriamo.
Prendiamo in esame un’equazione del tipo
AX = B,
dove A Rn×n, non singolare, X Rn×m matrice incognita, B Rn×m matrice dei termini noti.
Questo problema si riconduce facilmente al caso in cui esiste un solo vettore incognito e un solo
vettore dei termini noti. Infatti se X = (x1 x2 …xm), e B = (b1 b2 …bm), si ha AX = (Ax1 Ax2 …Axm) = (b1
b2 …bm). Questo perché esiste una classica interpretazione del prodotto tra due matrici A e X Il
prodotto consiste nel moltiplicare nell’ordine A per ogni vettore colonna di X e questi vettori,
nell’ordine formano le colonne di AX. Quindi, eguagliando colonna per colonna si trae Axi = bi, i
=1, …m. Si tratta di risolvere m sistemi lineari tutti con la stessa matrice dei coefficienti. Vedremo
che questi sistemi si possono risolvere uno alla volta dopo ave opportunamente “ manipolato” la
matrice dei coefficienti.
Un caso importante
Supponiamo che m = n e B = I. Osserviamo che la matrice identità I, enfatizzando le
colonne, è I = (e1 e2 …en). Allora
AX = I
equivale a risolvere gli n sistemi lineari
Axi = ei, i =1, …n.
L’unica matrice X soluzione di questa equazione matriciale è la matrice inversa di A, denotata con
X = A-1. Essa è dunque tale che
AA-1 = I
Se moltiplichiamo a destra per A si ha
A(A-1A) = IA = A
Posto Y = A-1 A, l’equazione AY = A, ammette la sola soluzione Y = I, da cui l’inversa di A verifica
anche A-1 A = I.
Se esiste una matrice A-1 che verifica una o l’altra delle due proprietà ( sono equivalenti) A si dice
invertibile. Allora una matrice è invertibile se e solo se è non singolare.
Soluzioni formali
Con la nozione di inversa possiamo formalmente scrivere la soluzione del sistema
AX = B
moltiplicando a sinistra per A-1
A-1AX = A-1B
Ottenendo, ( I X = X)
X = A-1B
Nel caso particolare di un solo vettore colonna di termini noti
x = A-1b.
Analogamente una equazione matriciale del tipo
XA = B
dove A𝜖 Rn×n non singolare, X 𝜖Rm×n matrice incognita, B 𝜖 Rm×n matrice dei termini noti si
risolve formalmente moltiplicando a destra per A-1,
XAA-1 = BA-1,
ottenendo ( X I = X)
X = BA-1.
OSSERVAZIONE
L’equazione matriciale XA = B si trasforma nella equivalente equazione matriciale
ATXT = BT.
Usando la regola del trasposto del prodotto tra matrici.
Allora, se
XT = (x
T1 x
T2 …x
Tm)
BT = (b
T1 b
T2 …b
Tm)
si tratta ancora una volta di risolvere m sistemi lineari
AT xTi = bT
i, i =1, …m
Questa volta il simbolo trasposto sui vettori colonna è stato usato abusivamente. Si vuole dire
semplicemente che sono i vettori colonna di XT e BT. Una volta calcolato XT, mediante
trasposizione si trova X.
Proprietà dell’inversa.
Se A, B , C, D, …..sono matrici quadrate di ordine n tutte non singolari tale è anche il prodotto
ABCD…..risultando
(ABCD….)-1 = …..D-1C-1B-1 A-1.
Poi , (A-1)-1 = A, ma, anche
(A-1)T = (AT)-1,
che giustifica la notazione A-T.
Regola di Binet
Risulta anche
det(ABCD….) = detA × detB × deC…..
Ricordiamo, anche, che det I = 1, e, in conseguenza,da ,detA-1
detA = detI = 1, si trae
detA-1 = 1/detA
Indice di condizionamento.
Riprendiamo in esame il sistema lineare
Ax = b
dove A è una matrice quadrata di ordine n. Se A è non singolare il sistema ammette la sola
soluzione x = A-1b. A è non singolare se e se detA ≠ 0. Dal punto di vista computazionale, diventa
importante capire cosa accade se A è “ quasi singolare”. Occorre in qualche modo introdurre una
misura che indichi quanto A sia distante dalla classe delle matrici singolari. Approssimativamente
sembrerebbe che A è “quasi singolare” se il suo determinante è “molto vicino a zero”. Non è
esattamente così. Scelta una norma definiamo, in corrispondenza, l’indice di condizionamento
come
κ(A) = || A || || A-1 ||
a patto che A sia non singolare. Ora, sappiamo che I = A A-1 e che, per ogni norma, || I || ≥ 1. (In
quelle indotte si ha || I || = 1). Per la quarta proprietà delle norme si ha
1 ≤ || I || = || AA-1|| ≤ ||A || || A-1 || = κ(A)
Cioè l’indice di condizionamento è sempre maggiore o uguale a uno. Ora moltiplichiamo A per
uno scalare arbitrariamente grande o piccolo c non nullo, e calcoliamo l’indice di
condizionamento di cA. Risulta (cA)-1 = A-1/ c, da cui
κ(cA) = || cA || || (cA)-1 || = | c| || A|| || A-1 || /|c| = || A || || A-1 ||= κ(A)
Segue che l’indice di condizionamento è invariante per prodotto con uno scalare. In sostanza non
varia se la matrice viene scalata o verso il basso o verso l’alto. L’indice di condizionamento detto
anc numero di condizione misura, nella norma scelta, il discostarsi o meno di A dalla classe matrici
singolari. A seconda della norma scelta si hanno le scritture κ 2(A), κ1(A), κ (A), ecc.Una matrice
in cui κ (A) = 1, ( il caso ideale ) è ottimamente condizionata , ma è ottimamente condizionata
anche se κ (A) = 7 o κ (A) = 10, e si ha un molto soddisfacente condizionamento se κ (A) = 100 o κ
(A) = 200. Probabilmente è da ritenersi ben condizionata una matrice con κ (A) =1000,o anche κ
(A) = 10000 ( dipende dai contesti). Una matrice comincia a diventare mal condizionata quando il
numero di condizione comincia ad essere elevate potenze di 10 come 107, 1010, 1020, ecc. In questi
casi la matrice A è molto vicina ad una matrice singolare. Per meglio comprendere riprendiamo
la regola di Laplace per l’inversa di A
1,...nj1,...n,i,detA
detA1)(a
ij
ji
1
ij
dove Aij si ottiene da A eliminando riga i e colonna j. Si vede che gli elementi dell’inversa
diventano tanto più grandi quanto più piccolo è il determinante di A. Ma, se questi elementi sono
molto grandi, molto grande diventa anche || A -1 ||. Sembrerebbe che il numero di condizione
diventa anch’esso molto grande. Ciò è soloo parzialmente vero, perché se ogni aij è molto piccolo,
molto piccola è anche ||A || e il prodotto || A || || A -1 || potrebbe rimanere in limiti accettabili.
Esempi.
Sia
21
11A
Si ha det A = 1, e
11
121A
Poi || A ||1 = 3, || A-1 ||1 = 3, e
κ1(A) = 9
La matrice A è molto ben condizionata. Lo stesso indice di condizionamento si ha in norma
infinito, mentre κ2(A) = 6.85… Ora moltiplichiamo A per c = 10-15, ottenendo
1515
1515
110210
1010A
e detA1 = 10-30 !! Un determinante così basso farebbe pensare ad una matrice mal condizionata,
anche perché la sua inversa diventa
1515
1515
1
11010
10102A
Questa inversa ha elementi grandi e si ha
151
1 103|||| A
Fortunatamente A1 ha elementi molto piccoli è e si ha ||A1|| = 9. Questo implica κ1(A1) = 9. Cioè
Anche A1 è una matrice ben condizionata !! . D’altra parte non c’è da meravigliarsi dal momento
che abbiamo già provato che l’indice di condizionamento è invariante per prodotti con scalari. Un
po’ grossolanamente esprimendoci possiamo dire che una matrice A è mal condizionata se ha un
determinante molto piccolo una volta che la si normalizzi in modo che l’elemento di modulo più
alto sia di qualche unità ( per esempio uno). Normalizzare vuol dire dividere A per una costante
opportuna.
Consideriamo ora la matrice
151011
11A
Si ha detA = 10-15. Questa volta il determinante piccolo è una spia di mal condizionamento. Infatti
1515
1515
1
1010
10101A
Poi, || A || = 2 + 10-15 > 2, ||A-1|| = 2 × 1015 +1 > 2 × 1015, da cui κ (A) > 4 × 1015 e la
matrice è mal condizionata.
La sensitività di un sistema lineare.
Per meglio comprendere il significato di indice di condizionamento riprendiamo il sistema
lineare
Ax = b
con A matrice quadrata non singolare di ordine n.
Proviamo a perturbare gli elementi di A e b in modo che A resti non singolare. Queste piccole
perturbazioni siano espresse dalla matrice δA e dal vettore δb. La soluzione x = A-1b del sistema
originario, viene a modificarsi in x = x + δx, in quanto soluzione di
(A + δA)(x + δx) = b + δb.
Dobbiamo fare riferimento alle perturbazioni relative espresse, sceltauna norma su vettori e
matrici, da ||||
|||| e
||||
||||
b
b
A
A.Usiamo lo stesso pedice per la norma su vettori e matrici immaginando
che una sia indotta dall’altra ( per esempio la norma infinito o la norma uno). O se non si tratta di
norme indotte pensiamo che siano norme compatibili, vale a dire ||Bu|| ≤ ||B|| ||u||, per ogni
matrice B e vettore u , come è il caso tra norma due su vettori e norma di Frobenius su matrici.
Allora, si dimostra che
||||
||||
||||
||||)(
||||
||||
b
b
x
x
A
AAc
dove c è una costante maggiore di uno , ma che è di qualche unità. Per esempio, se
2
1
||||
||||)(
A
AA
risulta c = 2 e per più piccoli valori di ||||
||||
A
Atale costante si attesta intorno a uno. La
disuguaglianza che mette in relazione la perturbazione sulla soluzione x e le perturbazioni sulla
matrice A e il vettore termine noto b, è stretta, vale a dire che minore o uguale può essere
sostituito con “ quasi uguale” e, per qualche matrice, direttamente col segno uguale. Allora
piccole perturbazioni relative ||||
|||| e
||||
||||
b
b
A
A vengono amplificate dal numero di condizione κ (A)
per ottenere la perturbazione relativa sul vettore soluzione ||||
||||
x
x. Se questo numero di
condizione è molto grande si hanno grosse perturbazioni relative sul vettore soluzione. Vale a dire
chse la matrice è “quasi singolare” perturbando di poco gli elementi di A b si hanno grosse
perturbazioni relative sulle componenti di x.
Conseguenze computazionali.
Supponiamo che nel risolvere il sistema
Ax = b
e l’unico errore commesso sia quello di trasformare gli elementi di A e b in numeri floating
dell’aritmetica finita usata. Cioè, supponiamo idealmente, di non commettere altri errori, di non
dovere più usare la funzione floating nelle operazioni aritmetiche. Allora se u è l’unità di
arrotondamento dell’aritmetica, si ha che in luogo della matrice A usiamo la matrice di coefficienti
fl(aij) = aij( 1 + εij ), = aij + aij εij, | εij | < u. Posto δA = fl(A) –A, si ha che | δA| < u | A |. Se usiamo
norme tali che la norma del valore assoluto di una matrice è uguale alla norma della matrice ( per
esempio norma uno e norma infinito), e valga la proprietà di monotonia , cioè O ≤ X ≤ Y, implica
||X || ≤ || Y ||, ( per esempio norma uno e norma infinito) , si trae || δA|| < u||A||. In modo
analogo se si usa fl(b), si trova, posto δb=fl(b) – b, || δ b || < u ||b|| . In definitiva
.||||
||||
||||
||||u
b
bu
,
A
A
In conseguenza se c = 2
ux
x)(4
||||
||||A
Leggiamo questa relazione ( il 4 non ha nessuna importanza) come
))((||||
||||u
x
xAO
Cioè l’ordine di grandezza dell’errore relativo sulla soluzione è il prodotto tra l’unità di
arrotondamento e il numero di condizione. Questo è detto errore minimo perché non si può
prescindere da questo errore qualunque algoritmo si usi per risolvere il sistema. L’errore reale
commesso è sempre più grande di tale errore. Per esempio in ambiente Matlab in cui u = eps =
O(10-16), ogni qualvolta si risolve un sistema Ax = b, si commette un errore che è sicuramente
maggiore di O(κ (A) × 10-16). Se κ (A) = 1013, l’errore minimo è 10-3, vale a dire al più tre cifre
decimali significative corrette!! ( rispetto alle 16 di partenza). Se κ (A) = 1017, non c’è alcuna
speranza di approssimare decentemente la soluzione, ( a meno di precondizionamenti), ecc.
I metodi diretti
I metodi di risoluzione che ci accingiamo ad illustrare sono sostanzialmente due. Il primo è
il metodo di fattorizzazione LU, e, il secondo e il metodo di fattorizzazione LU con la strategia del
pivoting. La parola fattorizzazione significa che si riesce a scrivere la matrice dei coefficienti come
prodotto di due matrici “ semplici”, esattamente L e U, cioè tali che i sistemi lineari ad esse
associate siano di “facile” risoluzione. Storicamente l’idea è quella, probabilmente familiare, di
triangolarizzare il sistema di equazioni in modo da calcolare le incognite una alla volta. Questo
procedimento è stato nel tempo rivisitato in termini di fattorizzazione per renderlo più adatto ai
problemi reali. Questi due metodi forniscono algoritmi ormai ben sviluppati in ottimi codici di
calcolo, ma non sono in grado di risolvere sistemi lineari di grosse dimensioni che intervengono ,
per esempio, nella discretizzazione di equazioni alle derivate parziali. Per bene introdurre le
metodologie, cominciamo con un esempio. Sia da risolvere il sistema
x1 + x2 + 3x3 = -5
2x1 –2x2 + x3 = 7
4x1+ 2x2 + 8x3 = -10
Approfittiamo del fatto che il coefficiente di x1 nella prima equazione è non nullo, Esso vale 1 e
lo battezziamo primo pivot. La prima equazione è l’equazione pivotale, e la prima riga della
matrice dei coefficienti è la riga pivotale. Per eliminare l’incognita x2 dalla seconda equazione,
usiamo il moltiplicatore l21 = a21/a11 = a21/ pivot = 2/1 = 2. Quindi sottraiamo dalla seconda
equazione la prima equazione premoltiplicata per il moltiplicatore l21. Questo vuol dire che prima
di fare la differenza tra coefficienti corrispondenti, quelli della equazione pivotale devono essere
moltiplicati per l21 = 2. Poi, eliminiamo l’incognita x3 dalla terza equazione, usando il moltiplicatore
l31 = a31/a11 = a31/pivot = 4/1 = 1, e sottraendo dalla terza equazione l’equazione pivotale
premoltiplicata per il moltiplicatore l31. In questo modo, termina il primo passo, conducendo al
sistema equivalente
x1 + x2 + 3x3 = -5
-4x2 –5x3 = 17
-2x2 –4x3 = 10
I moltiplicatori li1, i =2, i =3, hanno come secondo indice un numero che è quello dell’equazione
pivotale o della riga pivotale o del passo, e primo indice quello delle righe successive ( o
equazioni successive). Dimentichiamo, quindi, la prima equazione soffermandoci sul sistema di
due equazioni in due incognita individuata dalla seconda e terza equazione. Approfittiamo del
fatto che il coefficiente dell’incognita x2 nella seconda equazione è non nullo. Esso vale -4 e
diventa il secondo pivot . Adesso la seconda equazione è l’equazione pivotale e la seconda riga
della ( nuova) matrice dei coefficienti è la riga pivotale. Per eliminare l’incognita x3 dalla terza
equazione , costruiamo il moltiplicatore a32/a22 = a32/pivot = -2/-4 = ½. (l’elemento a22 non è
quello della matrice originaria). Sottraiamo, quindi, dalla terza equazione l’equazione pivotale ,
premoltiplicata per il moltiplicatore l32 . Questo vuol dire che prima di fare la differenza tra
elementi corrispondenti delle due equazioni, quelli della equazione pivotale vanno moltiplicati per
l32. In questo modo termina il secondo passo. Il moltiplicatore l32 ha come secondo indice un
numero che è quello della equazione pivotale o riga pivotale o del passo, mentre il primo indice è
quello della riga successiva o equazione successiva. Concluso questo secondo passo, il sistema si
trasforma in
x1 + x2 + 3x3 = -5
-4x2 –5x3 = 17
(-3/2)x3 = 3/2.
Il nuovo sistema di equazioni lineari è del tutto equivalente al primo, salvo che è più facilmente
risolubile. Dalla terza equazione, possiamo calcolare x3 = -1, e sostituendo nella seconda, possiamo
calcolare x2 = -3. Infine sostituendo i valori calcolati di x1 e x2, nella prima equazione, siamo in
grado di calcolare anche x1 = 1. In definitiva il vettore
x = (1 –3 –1)T
è il vettore soluzione del sistema proposto e questo si controlla a posteriori, sostituendo i valori
della tre incognite nel sistema originario. Da notare che il numero dei passi impiegati è stato 3 – 1
= 2, cioè uno in meno della dimensione. Supponiamo, ora, si voglia risolvere un nuovo sistema di
equazioni lineari con la stessa matrice dei coefficienti, ma con un altro vettore termine noto. Per
esempio, il nuovo sistema sia
x1 + x2 + 3x3 = 5
2x1 –2x2 + x3 = 1
4x1 + 2x2 + 8x3 = 14
Se non conservassimo i moltiplicatori precedentemente creati, dovremmo ripetere le stesse
operazioni per trasformare il sistema di equazioni nella forma semplice prima costruita.
Naturalmente, la stessa situazione si verrebbe a creare con un terzo vettore termine noto, un quarto e
così via. Se i termini noti si conoscessero prima, basterebbe allora mettere i vettori termini noti uno
dopo l’altro a formare una matrice di termini noti. Quindi le operazioni precedenti si potrebbero
estendere a tutti i termini noti e risolvere uno alla volta i sistemi lineari finali prendendo un vettore
alla volta. Ma questo non è sempre il caso. Per esempio, può capitare che il nuovo vettore termine
noto sia il vettore incognito precedente o una funzione del vettore incognito precedente. Allora, una
volta deciso di bene utilizzare i moltiplicatori, organizziamoli nella semplice matrice , con elementi
diagonali tutti uguali a uno.
11/24
012
001
L
1
01
001
3231
21
ll
l
I moltiplicatori sono stati posizionati rispettando i loro indici. Poi, dal sistema trasformato nella
forma finale, prendiamo la matrice dei coefficienti( che è triangolare superiore ).
3/200
540
311
u00
uu0
uuu
U
33
2322
131211
Infine consideriamo la matrice dei coefficienti del sistema originario
824
122
311
aaa
aaa
aaa
A
333231
232221
131211
Ci accorgiamo, allora, che risulta
A = LU
Una parte delle operazioni precedentemente eseguite, ha formato allora la fattorizzazione della
matrice dei coefficienti in una matrice L, che possiamo chiamare triangolare inferiore o triangolare
bassaa, e in una matrice U che possiamo chiamare triangolare superiore o triangolare alta. Da
notare che L ha anche elementi diagonali pari a uno. Le operazioni precedentemente descritte si
possono non eseguire quando riguardano i coefficienti del termine noto.
Infatti, scriviamo il sistema nella forma
LUx = b.
Per la proprietà associativa del prodotto tra matrici si può scrivere
L(Ux) = b.
Poniamo allora
Ux = y.
Il vettore y = (y1, y2, y3)T è soluzione del sistema
Ly = b
In forma esplicita, tale sistema si scrive
y1 = -5
2y1 + y2 = 7
4y1 + (1/2)y2 + y3 = -10
La prima equazione dà subito y1 = -5, e, sostituendo nella seconda, si trova y2 = 17; infine, con i
valori trovati di y1 e y2,
dalla terza equazione si trae y3 = 3/2. Ora, il vettore soluzione x è tale che
Ux = y,
cioè esplicitamente
x1 + x2 + 3x3 = -5
-4x2 –5x3 = 17
(-3/2)x3 =3/2.
Questa è proprio la forma finale del sistema precedentemente trovata, con soluzione x3 = -1, x2 =
3, x1 = 1. In sostanza la soluzione è stata trovata risolvendo due sistemi triangolari. Il primo con
matrice triangolare bassa e il secondo con matrice triangolare alta. Si dice che il primo sistema è
stato risolto con il metodo di sostituzione in avanti e il secondo con il metodo di sostituzione
all’indietro. Ma, se cambiamo il vettore termine noto, per esempio, b = (5 1 14)T, il sistema con
matrice triangolare bassa Ly = b, diventa
y1 = 5
2y1 + y2 = 1
4y1 + (1/2)y2 + y3 = 14
Allora, y1 = 5, e sostituendo nella seconda equazio y2 = -9, e con questi valori nella terza, si trae
infine y3 = -3/2. Il sistema con matrice triangolare alta Ux = y, diventa invece
x1 + x2 + 3x3 = 5
-4x2 –5x3 = -9
(-3/2)x3 = -3/2.
Dalla terza equazione si trova x3 = 1, e sostituendo nella seconda, si trova x2 = 1, e mettendo
entrambi questi valori nella prima, si ottiene ancora x1 = 1.
Il vettore soluzione è ora x = (1 1 1)T. Questo vettore è stato ottenuto riutilizzando gli elementi L e
U della fattorizzazione. Con un terzo, quarto vettore, ecc., come termine noto, si ottiene allo stesso
modo il vettore soluzione.
In sostanza conviene preliminarmente eseguire la fattorizzazione della matrice A e
successivamente, per ogni vettore termine noto trovare il corrispondente vettore soluzione. La
procedura che ha portato alla triangolarizzazione del sistema originario, è un esempio d riduzione
di Gauss, ( o riduzione gaussiana), mentre la interpretazione in termini di prodotto di due
matrici è un esempio di fattorizzazione di Gauss ( o fattorizzazione gaussiana). In onore del
grande matematico.
Questo algoritmo si estende a matrici di tutte le dimensioni come vedremo in seguito.
Osserviamo, che, questo problema è stato risolto, perché entrambi i pivot incontrati erano non
nulli e permettevano di formare i moltiplicatori. Nella generalizzazione che andremo a fare per
matrici qualsiasi, questo, in generale, non è vero. Il fatto che si può incontrare un pivot nullo è il vero limite delle procedure gaussiane usate senza
particolari accorgimenti.