Calcolo numerico - Aappunti

P

det(A) = fisso una i = a (−1) det(A ) =

ij ij

j=1

i+1 i+2 i+n

a (−1) det(A ) + a (−1) det(A ) + . . . + a (−1) det(A )

i1 i1 i2 i2 in in

dove per A si intende la sottomatrice ottenuta da A eliminando la riga i e la colonna 1.

i1

Analogamente la regola vale se fisso una j e quindi procedo per colonne.

Si rammenta che la regola di Laplace non è direttamente utile per risolvere il sistema, ma necessaria per la

regola di Kramer.

Regola di Kramer ↔ ⃗b

′

nxn n n

∈ ∈ ∈ ∋

Sia A , b , x A⃗x = allora risulta che:

R R R det(A )

∀i i

= 1, 2, 3, . . . , n x =

i det(A)

dove A è la matrice ottenuta cambiando la i-esima colonna con quella dei termini noti (b).

i

Regola di Kramer costo computazionale

Per ogni i devo calcolare det(A ) e cioè lo calcolo n volte in più devo calcolare det(A) una volta. In totale:

i

(n+1)*C(det). Con dei piccoli esempi mi rendo conto che il metodo di Kramer è impraticabile poiché per

157

esempio per n = 100 ottengo C = 10 .

Sistemi lineari triangolari

Supponiamo di avere:      

a b c x g

1 =

0 d e x h

     

2

     

0 0 f x i

3

notiamo che A è una matrice triangolare superiore.

    

 x g

a 0 0 1 = h

b c 0 x

     

2

    

 d e f x i

3

notiamo che A è una matrice triangolare inferiore. 10

Tecnica di sostituzione all’indietro o all’avanti

Se ho una matrice triangolare superiore procedo con la sostituzione all’indietro:

• b ;

x = n

n a nn n

P

−

b a x

• i ij j

j=i+1 ∀i −

x = ; = n 1, . . . , 1.

i a

ii

Se ho una matrice triangolare inferiore procedo con la sostituzione in avanti:

• b

x = ;

1

1 a

11 i−1

P

−

b a x

• i ij j

j=1

x = ; i = 2, . . . , n.

i a ii

Costo computazionale sistema triangolare

n 2

P −

C(ST) = (2i 1) = n .

i=1

Basta considerare che 1 (divisione) + 1 (sottrazione) + (i-2) (sommatoria) + (i-1) prodotto è il numero di

operazioni per ogni i (cioè 2i-1).

Metodo di eliminazione di Gauss

Col metodo di eliminazione di Gauss passiamo da un sistema lineare quadrato ad uno triangolare.

nxm n n

∈ ∈ ∈

Sia A , x , b se Ax=b:

R R R ↔ ↔ ↔ ↔ ↔

B⃗b

↔ ↔ ↔

Ax = b B A⃗x = . . . T ⃗x = ⃗c

Le seguenti operazioni non modificano il sistema:

1. scambiare le righe del sistema;

2. moltiplicare e dividere per la stessa quantità entrambi i membri di una equazione;

3. somma membro a membro due equazioni.

Importante ricordare che nel MEG si usano i moltiplicatori di Gauss indicati con l mentre l’elemento

ij

che si vuole eliminare si chiama elemento pivotale.

Osservazioni MEG

• il MEG è iterativo;

• ha n-1 iterazioni;

• ∀k −

= 0, 1, . . . , n 1 : rimane costante la riga k-1, modifico la colonna k (sotto diagonale).

k k ∀k −

Il metodo di eliminazione di Gauss crea una successione di coppie A e b tali che = 0, 1, . . . , n 1 si

k k

abbia: A x = b che sono tutti sistemi equivalenti.

Gli steps sono:

• a

−

l = ;

ik

ik a kk 11

• l (rigak);

in

• riga i + l (rigak)

in

che equivalgono a: k k

n k k

a a

k k

P − −

a a = b b .

ik ik

ij kj i k

k k

j=i+1 a a

kk kk

Costo computazionale MEG

Per il calcolo del costo computazionale del metodo di eliminazione di Gauss dobbiamo considerare il costo

computazionale per passare da un sistema quadrato ad uno triangolare + il costo computazionale per risolvere

2

il sistema triangolare; quest’ultimo sappiamo essere pari ad n per cui procediamo col calcolare il primo.

Considero: k k

1. numero di operazioni per modificare 1 elemento di A e b

k k+1 k k+1

→ →

2. numero di operazioni totali per passare da A A e da b b ;

0 n−1 0 n−1

→ →

3. numero di operazioni per passare da A A e b b .

Con un paio di considerazioni si arriva a: n−1

0 n−1 0 n−1 k k+1 k k+1

P

→ → → →

C(A A , b b ) = C(A A , b b )

k=1

2

2 n 76

3 −

si ricava poi che quest’espressione vale: n + n a questo punto sappiamo il costo computazionale per

3 2

passare da una matrice quadrata ad una triangolare e quindi rimettendo assieme i pezzi si ha:

23 32 76

3 2 −

C(M EG + T S) = n + n n

3

che si approssima a n .

Osservo che: C(MEG+TS) << C(KRAMER)

2 3

n << n!.

3

Strategie di pivoting

Servono ad evitare di avere una matrice singolare.

• Pivoting parziale:

krk kik

|a | | ≤ ≤

1. determinare = max|a con k i n;

2. scambio le equazioni r e k. 12

Nella zona rosa cerco il più grande in modulo e sarà a dopodiché scambio riga r con riga k.

rk

• Pivoting totale:

krs kij

|a | |a |,

Sia ora a , = max un qualsiasi elemento nella zona rosa scambio le righe k ed r e scambio

rs k≤i,j≤n

le colonne k ed s per ottenere un pivoting totale.

Fattorizzazione LU

Supponiamo un problema che ad un determinato step richieda la risoluzione con il metodo di Gauss del

sistema lineare Ax = b; la matrice diviene triangolare superiore e quindi viene richiesto di risolvere il sistema

n n

triangolare A x = b se si ipotizza che nello stesso problema dopo un certo tempo sia necessario risolvere il

sistema Ax = c che ha la stessa matrice A del precedente, ma diverso vettore dei termini noti (c al posto di

b) ci rendiamo conto che non possiamo usare i conti già fatti e per cui siamo costretti a dover riutilizzare

n n

Gauss per poi risolvere il sistema triangolare A x = c . L’argomento di questo paragrafo ci presenterà un

metodo per evitare di dover rifare tutti i calcoli.

Si stabilisce l’esistenza di una matrice L, triangolare inferiore e con elementi della diagonale pari a 1, e di

una matrice U, triangolare superiore, tali che: A = LU

ciò implica che: det(A) = det(L)det(U ), ma det(L) = 1 per cui det(A) = det(U ).

A questo punto possiamo dire che il sistema di partenza Ax = b si può scrivere come LUx = b. Per questo

↔ ⃗ ⃗

⃗b

′ ′

motivo ora mi ritrovo a dover risolvere il sistema L x = e dopo aver ricavato x vado a risolvere il sistema

↔ ⃗

′

U⃗x = x .

Il costo computazionale per un solo vettore dei termini noti (b) è:

23 3

≃ n .

C(LU)

Ricorda che applicando il MEG al sistema (Ax = b) vuol dire ottenere una matrice triangolare superiore

n−1 n−1 n−1

(A ) e cioè il nuovo sistema A x = b che è un sistema triangolare; mentre quando applico la

fattorizzazione LU al sistema (Ax = b) vuol dire ottenere due sistemi: il sistema Ly = b ed il sistema Ux =

y e cioè due sistemi triangolari.

Proprietà n

• P

→

A = LU a = l u , ma se L ed U sono triangolari il prodotto riga per colonna è pari a:

ij ik kj

k=1

min(i,j)

P → ∀i ∈ ∀j ∈ →

* l u l = 1 [1, n] allora : a = u [1, n] prima riga di A = prima riga

ik kj ii 1j 1j

k=1

di U. 13

min(i,j) i i−1

• P P P

≤

Se i j : a = l u = l u = l u + l u = l valgono 1 = u +

ij ik kj ik kj ik kj ii ij ii ij

k=1 k=1

k=1 i−1

i−1 P

P −

→ u = a l u .

l u ij ij ik kj

ik kj k=1

k=1 j−1

P

−

a l u

j j−1

• ij ik kj

P P →

Se i > j : a = l u = l u + l u .

l = k=1

ij ik kj ik kj ij jj ij

k=1 k=1 u jj

∈ −

Osservo che per calcolare l devo conoscere tutti gli u con k [1, i 1], mentre per calcolare u devo

ij kj ij

∈ −

conoscere tutti gli l con k [1, j 1].

ik

Spiegazione * :

Tecnica di Crout (vado avanti per righe)

Steps: ∧

1. l = 1 u = a per i = 1;

11 1j 1j

2. calcolo l , l è 1, calcolo u per j∈ [2, n];

21 22 2j

∈ − ∈

k. calcolo l per j [1, k 1], l è 1, calcolo u per j [k, n];

kj kk kj

∈ −

n. calcolo l per j [1, n 1], l è 1, calcolo u .

nj nn nn

Rappresentazione:

i numeri indicano gli step.

(Procedo calcolando la prima riga di U e la seconda riga di L e cosı̀ via).

14

Tecnica di Doolittle (mix righe colonne)

Per le righe ragioniamo con U e per le colonne con L.

Steps:

1. calcolo prima riga a = u ;

ij ij

2. calcolo prima colonna l ;

i1

3. calcolo la seconda riga con gli elementi u ;

2j

4. calcolo seconda colonna l ;

i2

k. calcolo la riga k-esima u ;

kj

k+1. calcolo la colonna k-esima l ;

ik

n. calcolo elemento finale u .

nn

Rappresentazione:

(Procedo calcolando alternamente la prima riga di U e la prima colonna di L e cosı̀ via).

Osservazioni fattorizzazione LU

• no pivoting in LU;

• Crout equivale a Doolittle.

Equivalenza tra MEG ed LU

Con il metodo di Gauss si va a costruire una successione di sistemi lineari equivalenti (per l’esattezza n-1

sistemi). ′ ′ ′′ ′′

0 0 n−1

→ → → →

Analizziamo il processo: dato il problema Ax = b A x = b A x = b A x = b A x =

n−1 n−1 n−2 n−3 0 n−1 n−2 n−2

→ →=

b det(A ) = det(A ) = det(A ) = . . . = det(A ) = det(A) det(A ) = det(L A ) =

n−2 n−2

T h. Binet = det(L )det(A ).

Abbiamo: 15

dove gli elementi l sono i moltiplicatori di Gauss che calcolo nel primo step di MEG.

La matrice L della fattorizzazione LU è composta da tutti i moltiplicatori di Gauss con segno opposto.

n−1

La matrice U della fattorizzazione LU è proprio A di MEG.

Condizionamento di un sistema lineare

Ci si è resi conto nel primo capitolo che a causa dell’aritmetica finita dell’elaboratore i risultati si possono

considerare inaffidabili, proprio per questo motivo alcune volte il calcolo potrebbe fornire dei risultati molto

lontani da quelli reali.

Si consideri il seguente sistema lineare:

Esempio x + y =2

1000x + 1001y = 2001

la soluzione è x = y = 1. Si pensi ora di perturbare dell’1% il coefficiente x della prima equazione, il sistema

diviene: (1 + 0.01)x + y = 2

1000x + 1001y = 2001

benché si possa immaginare che la soluzione del sistema non sia molto diversa da quella del precedente in

1901

1

− e y = .

questo caso essa risulta essere x = 9 900

Un sistema lineare per cui