Algebra lineare numerica

Alvise Sommariva

Università degli Studi di Padova Dipartimento di Matematica Pura e Applicata

30 novembre 2011

Introduzione

Sia A ∈ ℝ^{n×n} una matrice a coeff. reali, b ∈ ℝ^{n} un vettore colonna e supponiamo di dover calcolare un vettore colonna x^{∗} ∈ ℝ^{n} cosicchè A · x^{∗} = b.

Come è noto questo problema ha soluzione unica x^{∗} se e solo se det (A) ≠ 0 (matrice non singolare). Ci porremo di seguito in queste ipotesi.

Introdurremo due famiglie di metodi, quelli diretti che calcolano la soluzione x^{∗} in un numero finito di operazioni e quelli iterativi per cui la soluzione x è ottenuta come limite di un processo iterativo.

Eliminazione gaussiana: esempio.

Risolvere il sistema lineare

1 · x_{1} + 2 · x_{2} + 1 · x_{3} + 4 · x_{4} = 13 2 · x_{1} + 0 · x_{2} + 4 · x_{3} + 3 · x_{4} = 28 4 · x_{1} + 2 · x_{2} + 2 · x_{3} + 1 · x_{4} = 20 −3 · x_{1} + 1 · x_{2} + 3 · x_{3} + 2 · x_{4} = 6

utilizzando esclusivamente le seguenti regole, che modificano il sistema lineare da risolvere ma non le sue soluzioni:

▶ Scambi di righe; ▶ Moltiplicazioni di righe per una costante; ▶ Una riga può essere sostituita con la somma della stessa e un

multiplo non nullo di un’altra riga.

Eliminazione gaussiana: esempio.

Posti

A = ( \begin{matrix}1 & 2 & 1 & 4 \\ 2 & 0 & 4 & 3 \\ 4 & 2 & 2 & 1 \\ −3 & 1 & 3 & 2 \end{matrix} ) , b = ( \begin{matrix}13 \\ 28 \\ 20 \\ 6 \end{matrix} )

risolviamo il sistema lineare Ax = b con il metodo dell’eliminazione gaussiana. Consideriamo pure la matrice aumentata Aaug^{(1)} = (a_{j,k}^{(1)}):

Aaug^{(1)} = [ \begin{matrix}1 & 2 & 1 & 4 & 13 \\ 2 & 0 & 4 & 3 & 28 \\ 4 & 2 & 2 & 1 & 20 \\ −3 & 1 & 3 & 2 & 6 \end{matrix} ] .

Siano m_{k,1} = a_{k,1}^{(1)}/a_{1,1}^{(1)} i cosidetti moltiplicatori. L’elemento a_{1,1} è detto pivot.

\begin{matrix}m_{1,1} = 1 \\ m_{2,1} = 2 \\ m_{3,1} = 4 \\ m_{4,1} = −3 \end{matrix} [ \begin{matrix}1 & 2 & 1 & 4 & 13 \\ 2 & 0 & 4 & 3 & 28 \\ 4 & 2 & 2 & 1 & 20 \\ −3 & 1 & 3 & 2 & 6 \end{matrix} ] .

Per k = 2, 3, 4, moltiplichiamo la prima riga per −m_{k,1} e sommiamo il risultato alla k-sima riga. Abbiamo:

Aaug^{(2)} = [ \begin{matrix}1 & 2 & 1 & 4 & 13 \\ 0 & −4 & 2 & −5 & 2 \\ 0 & −6 & −2 & −15 & −32 \\ 0 & 7 & 6 & 14 & 45 \end{matrix} ] .

Siano m_{k,2} = a_{k,2}^{(2)}/a_{2,2}^{(2)} (k > 2) i moltiplicatori . L’elemento a_{2,2} è un pivot.

\begin{matrix}m_{2,2} = 1 \\ m_{3,2} = 1.5 \\ m_{4,2} = − − 1.75 \end{matrix} [ \begin{matrix}1 & 2 & 1 & 4 & 13 \\ 0 & −4 & 2 & −5 & 2 \\ 0 & −6 & −2 & −15 & −32 \\ 0 & 7 & 6 & 14 & 45 \end{matrix} ] .

Per k = 3, 4, moltiplichiamo la seconda riga per −m_{k,2} e sommiamo il risultato alla k-sima riga. Abbiamo:

Aaug^{(3)} = [ \begin{matrix}1 & 2 & 1 & 4 & 13 \\ 0 & −4 & 2 & −5 & 2 \\ 0 & 0 & −5 & −7.5 & −35 \\ 0 & 0 & 9.5 & 5.25 & 48.5 \end{matrix} ] .

Eliminazione gaussiana: esempio.

Siano m_{4,3} = a_{4,3}^{(3)}/a_{3,3}^{(3)} il moltiplicatore. L’elemento a_{3,3} è un pivot.

\begin{matrix}m_{3,3} = 1 \\ m_{4,3} = −1.9 \end{matrix} [ \begin{matrix}1 & 2 & 1 & 4 & 13 \\ 0 & −4 & 2 & −5 & 2 \\ 0 & −6 & −2 & −15 & −32 \\ 0 & 7 & 6 & 14 & 45 \end{matrix} ] .

Per k = 3, 4, moltiplichiamo la terza riga per −m_{k,3} e sommiamo il risultato alla k-sima riga. Abbiamo:

Aaug^{(4)} = [ \begin{matrix}1 & 2 & 1 & 4 & 13 \\ 0 & −4 & 2 & −5 & 2 \\ 0 & 0 & −5 & −7.5 & −35 \\ 0 & 0 & 0 & −9 & −18 \end{matrix} ] .

In virtù delle operazioni utilizzate, il sistema lineare

1 · x_{1} + 2 · x_{2} + 1 · x_{3} + 4 · x_{4} = 13
2 · x_{1} + 0 · x_{2} + 4 · x_{3} + 3 · x_{4} = 28
4 · x_{1} + 2 · x_{2} + 2 · x_{3} + 1 · x_{4} = 20
−3 · x_{1} + 1 · x_{2} + 3 · x_{3} + 2 · x_{4} = 6

ha quindi le stesse soluzioni del sistema lineare

1 · x_{1} + 2 · x_{2} + 1 · x_{3} + 4 · x_{4} = 13
0 · x_{1} + −4 · x_{2} + 2 · x_{3} + −5 · x_{4} = 2
0 · x_{1} + 0 · x_{2} − 5 · x_{3} − 7.5 · x_{4} = −35
0 · x_{1} + 0 · x_{2} + 0 · x_{3} − 9 · x_{4} = −18

che si risolve facilmente tramite sostituzione all’indietro (si calcola facilmente dall’ultima eqz. x_{4} = 2, si sostituisce tale valore alla penultima eqz. e si ottiene x_{3} = 4, si sostituiscono tali x_{3}, x_{4} alla seconda eqz. e si ottiene x_{2} = −1, , si sostituiscono tali x_{2}, x_{3}, x_{4} alla prima eqz. e si ottiene x_{1} = 3).

Alvise Sommariva  Algebra lineare numerica  8/ 25

Fattorizzazione LU

Il proposito è di rivedere questo algoritmo dal punto di vista matriciale.

La matrice iniziale

A = ( \begin{matrix}1 & 2 & 1 & 4 \\ 2 & 0 & 4 & 3 \\ 4 & 2 & 2 & 1 \\ −3 & 1 & 3 & 2 \end{matrix} )

si scrive come prodotto A = LU di due matrici quadrate L, U con

▶ L = (l_{i,j}) matrice triangolare inferiore, cioè l_{i,j} = 0 se i < j,

con la proprietà che l_{i,i} = 1. ▶ U = (u_{i,j}) matrice triangolare superiore, cioè l_{i,j} = 0 se i > j.

Introduciamo la routine fattLU:

f u n c t i o n B = fattLU ( A )
n=s i z e ( A , 1 ) ;
f o r k=1:n−1
A ( k+1:n , k )=A ( k+1:n , k ) /A ( k , k ) ;
f o r j=k+1:n ,
f o r s=k+1:n
A ( s , j )=A ( s , j )−A ( s , k ) ∗A ( k , j ) ;

end end

end

B=A ;

Esegue, se esiste, la fattorizzazione LU della matrice, immagazzinando opportunamente le matrici L, U nell’output B.

>> A=[1 2 1 4 ; 2 0 4 3 ; 4 2 2 1 ; −3 1 3 2 ] ;

>> B = fattLU ( A )
B =

1.0000  2.0000  1.0000  4.0000
2.0000  −4.0000  2.0000  −5.0000
4.0000  1.5000  −5.0000  −7.5000
−3.0000  −1.7500  −1.9000  −9.0000

>>

Si vede che il sistema lineare finale è immagazzinato nella parte diagonale superiore di B, mentre la parte diagonale inferiore contiene i moltiplicatori.

Fattorizzazione LU

Usando il comando triu che estrae la parte triangolare superiore della matrice:

>> U=t r i u ( B )

U =

1.0000  2.0000  1.0000  4.0000
0  −4.0000  2.0000  −5.0000
0  0  −5.0000  −7.5000
0  0  0  −9.0000
>> L=B−U+eye ( s i z e ( B ) )

L =

1.0000  0  0  0
2.0000  1.0000  0  0
4.0000  1.5000  1.0000  0
−3.0000  −1.7500  −1.9000  1.0000
>> norm ( A−L∗U )
ans =

0 >>

Alvise Sommariva  Algebra lineare numerica  12/ 25

Se A = LU allora, posto y := Ux abbiamo

Ly = LUx = Ax = b

Questo è un sistema triangolare inferiore che può essere agevolmente risolto con il metodo della sostituzione all’indietro. Una volta ottenuto y^{∗}, osservato che Ux = y^{∗} con un metodo simile, detto di sostituzione in avanti si ottiene la soluzione x^{∗}.

Si vede che se n è la dimensione della matrice A allora, in generale, servono O(n^{3}/3) per calcolare x^{∗}.

Alcuni problemi:

▶ Non è detto esista la fattorizzazione LU di una matrice A. Ad esempio, come si vede subito con facili conti, la matrice

A = ( \begin{matrix}0 & 1 \\ 1 & 0 \end{matrix} )

non è fattorizzabile LU. ▶ La presenza di moltiplicatori m_{i,j} grandi può causare errori

nella risoluzione del sistema lineare. ▶ Il calcolo dei moltiplicatori m_{s,k} = a_{s,k}^{(k)}/a_{k,k}^{(k)}, visto in

precedenza nella risoluzione di un particolare sist. lin. si interrompe se qualche denominatore a_{k,k}^{(k)} è nullo.

Matrice di permutazione

Una matrice P si dice di permutazione se ha per ogni riga e colonna una sola componente non nulla ed uguale a 1.

Alcune proprietà ▶ P^{−1} = P^{T}. ▶ Se P_{i,j} = 1 allora la i-sima riga di PA è uguale alla j-sima riga

di A.

>> P=[0 0 1 ; 1 0 0 ; 0 1 0 ]
P = 0  0  1
1  0  0
0  1  0
>> A=[1 2 3 ; 4 5 6 ; 7 8 9 ]
A =1  2  3
4  5  6
7  8  9
>> P∗A
ans = 7  8  9
1  2  3
4  5  6
Alvise Sommariva  Algebra lineare numerica  15/ 25

Cerchiamo l’elemento più grande in modulo della prima colonna di A

[ \begin{matrix}1 & 2 & 1 & 4 \\ 2 & 0 & 4 & 3 \\ 4 & 2 & 2 & 1 \\ −3 & 1 & 3 & 2 \end{matrix} ] .

e scambiamo la riga corrispondente (la prima riga viene scambiata con la terza). Otteniamo così

A^{(1)} = [ \begin{matrix}4 & 2 & 2 & 1 \\ 2 & 0 & 4 & 3 \\ 1 & 2 & 1 & 4 \\ −3 & 1 & 3 & 2 \end{matrix} ] .

Siano m_{k,1} = a_{k,1}^{(1)}/a_{1,1}^{(1)} i cosidetti moltiplicatori. L’elemento a_{1,1} è detto pivot. Per k = 2, 3, 4, moltiplichiamo la prima riga per −m_{k,1} e sommiamo il risultato alla k-sima riga. Abbiamo:

\begin{matrix}m_{1,1} = 1 \\ m_{2,1} = 0.5 \\ m_{3,1} = 0.25 \\ m_{4,1} = −0.75 \end{matrix} [ \begin{matrix}4 & 2 & 2 & 1 \\ 0 & −1 & 3 & 2.5 \\ 0 & 1.5 & 0.5 & 3, 75 \\ 0 & 2.5 & 4.5 & 2.75 \end{matrix} ] .

Osserviamo che se la matrice è non singolare, l’elemento di massimo modulo della prima colonna è non nullo (altrimenti avrei colonna di zeri!).

Partiamo da:

[ \begin{matrix}4 & 2 & 2 & 1 \\ 0 & −1 & 3 & 2.5 \\ 0 & 1.5 & 0.5 & 3, 75 \\ 0 & 2.5 & 4.5 & 2.75 \end{matrix} ] .

Osserviamo che se la matrice è non singolare, l’elemento di massimo modulo della seconda colonna con indice di riga ≥ 2 è non nullo (altrimenti avrei due colonne proporzionali e la matrice sarebbe singolare). Scambiamo la seconda riga con l’indice di riga per cui si ottiene tale massimo (quarta riga). Abbiamo così

A^{(2)} = [ \begin{matrix}4 & 2 & 2 & 1 \\ 0 & 2.5 & 4.5 & 2.75 \\ 0 & 1.5 & 0.5 & 3, 75 \\ 0 & −1 & 3 & 2.5 \end{matrix} ] .

I moltiplicatori m_{k,2} = a_{k,2}^{(2)}/a_{2,2}^{(2)} della matrice sono m_{3,2} = 1.5/2.5 = −0.6, m_{4,2} = −1/2.5 = −0.4.

\begin{matrix}m_{2,2} = 1 \\ m_{3,2} = −0.6 \\ m_{4,2} = −0.4 \end{matrix} [ \begin{matrix}4 & 2 & 2 & 1 \\ 0 & 2.5 & 4.5 & 2.75 \\ 0 & 1.5 & 0.5 & 3, 75 \\ 0 & −1 & 3 & 2.5 \end{matrix} ] .

Per k = 3, 4, moltiplichiamo la seconda riga per −m_{k,1} e sommiamo il risultato alla k-sima riga. Abbiamo:

[ \begin{matrix}4 & 2 & 2 & 1 \\ 0 & 2.5 & 4.5 & 2.75 \\ 0 & 0 & −2.2 & 2.1 \\ 0 & 0 & 4.8 & 3.6 \end{matrix} ] .

Con le stesse idee, per quanto riguarda la terza colonna, scambiamo la quarta riga con la terza e abbiamo

A^{(3)} = [ \begin{matrix}4 & 2 & 2 & 1 \\ 0 & 2.5 & 4.5 & 2.75 \\ 0 & 0 & 4.8 & 3.6 \\ 0 & 0 & −2.2 & 2.1 \end{matrix} ] .

Il moltiplicatore è m_{4,3} = −2.2/4.8 ≈ −0.4583. Per k = 4, moltiplichiamo la terza riga per −m_{4,3} e sommiamo il risultato alla quarta riga. Abbiamo la matrice triangolare superiore U = A^{(4)}:

\begin{matrix}m_{3,3} = 1 \\ m_{4,3} = −2.2/4.8 \end{matrix} [ \begin{matrix}4 & 2 & 2 & 1 \\ 0 & 2.5 & 4.5 & 2.75 \\ 0 & 0 & 4.8 & 3.6 \\ 0 & 0 & 0 & 3.75 \end{matrix} ] .

▶ Abbiamo scambiato la prima riga con la terza.

▶ Abbiamo scambiato la seconda riga con la quarta.

▶ Abbiamo scambiato la terza riga con la quarta. Quindi è come scambiare la prima riga con la terza, la seconda con la quarta, la terza con la seconda e la quarta con la prima. Così

P = [ \begin{matrix}0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \end{matrix} ] .

Per come abbiamo effettuato gli scambi, tutti i moltiplicatori sono in modulo m_{k,s} ≤ 1.

I moltiplicatori sono ▶ m_{1,1} = 1, m_{2,1} = 0.5, m_{3,1} = 0.25, m_{4,1} = −0.75 ▶ m_{2,2} = 1, m_{3,2} = 0.6, m_{4,2} = −1/2.5 = −0.4. ▶ m_{3,3} = 1, m_{4,3} = −2.2/4.8 ≈ −0.4583.

e costituiscono gli elementi della matrice L.

Fattorizzazione LU in Matlab/Octave

Con Matlab/Octave, verifichiamo i risultati ottenuti.

>> A=[1 2 1 4 ; 2 0 4 3 ; 4 2 2 1 ; −3 1 3 2 ] ;
>> [ L , U , P]= l u ( A )

L =

1.0000  0  0  0
−0.7500  1.0000  0  0
0.5000  −0.4000  1.0000  0
0.2500  0.6000  −0.4583  1.0000

U =

4.0000  2.0000  2.0000  1.0000
0  2.5000  4.5000  2.7500
0  0  4.8000  3.6000
0  0  0  3.7500

P =

0  0  1  0
0  0  0  1
0  1  0  0
1  0  0  0
>> norm ( P∗A−L∗U )
ans = 0
Alvise Sommariva  Algebra lineare numerica  23/ 25

Vediamo cosa fa Matlab/Octave. Si dimostra che ogni matrice A ∈ ℝ^{n×n} è fattorizzabile come

P · A = L · U

dove ▶ P è una matrice di permutazione. ▶ L = (l_{i,j}) è una matrice triangolare inferiore, cioè l_{i,j} = 0 se

i < j, con la proprietà che l_{i,i} = 1. ▶ U = (u_{i,j}) è una matrice triangolare superiore, cioè l_{i,j} = 0 se

i > j. Quindi se PA = LU e P^{−1} = P^{T} allora A = (P^{T} ∗ L)U Se PA = LU e Ax = b, allora LUx = PAx = Pb e quindi posto y = Ux, prima si risolve con la sostituzione all’indietro Ly = Pb e quindi, se y^{∗} è la soluzione, si risolve con la sostituzione in avanti Ux = y^{∗}.

>> A=[1 2 1 4 ; 2 0 4 3 ; 4 2 2 1 ; −3 1 3 2 ] ;
>> [ L , U , P]= l u ( A ) ;
>> b =[13; 2 8 ; 2 0 ; 6 ] ;
>> c=P∗b ;
>> y=L\c ;
>> x=U\y

x =

3.0000 −1.0000 4.0000 2.0000

>>