Appunti completi di Metodi numerici - Parte 2

T

inferiore e da per scontato che nella parte sopra della matrice ci sia la stessa cosa. La U sarà la L . Quindi

avremo:

A = L U

Che si scrive anche:

T

A = L L

In forma matriciale sarà: ?

∑ < >

= =

=@A

In forma semplificata sarà: 97

Quindi avremo la L e la U che sono: B B ⋯ B ⋯ B

B 0 0 0 0 0 ⎡ ⎤

⎡ ⎤

B 0 0 0 0

B 0 B ⋯ B ⋯ B

⎢ ⎥

⎢ ⎥

⋯ ⋯ ⋱ 0 0 0 ⎢ ⎥

0 0 ⋱ ⋯ ⋯ ⋯

⎢ ⎥ ⎢ ⎥

B B ⋯ B 0 0 ⋯ B

0 0 0 B

⎢ ⎥ ⎢ ⎥

L = U =

⋯ ⋯ ⋯ ⋯ 0 0

⎢ ⎥ 0 0 0 0 ⋱ ⋯

⎢ ⎥

B B ⋯ B ⋱ B

⎣ ⎦ 0 0 0 0 0 B

⎣ ⎦

Se ad esempio vado a moltiplicare la seconda riga di L per la quarta colonna di U, vediamo che sia per la riga

che per la colonna sto facendo delle operazioni inutili in cui ci sono gli zeri, allora vorrà dire che sia per la riga

e per la colonna scelta mi devo fermare appena incontro uno zero, infatti nella formula precedente nella

T

sommatoria abbiamo min (i,j). Siccome la U = L allora nella formula non avremo u ma l . Quindi la formula

kj jk

T

in rosso è la formula di A = L L in termini matriciali però implementata in maniera furba da non fare calcoli

inutili quando ci sono gli zeri.

Quindi noi dalla forma matriciale:

T

A = L L

Abbiamo noti gli elementi di A, e vogliamo andarci a ricavare i valori degli elementi di L, così di conseguenza

T

riusciamo a ricavarci anche quelli di L . Quindi prendo un elemento di A e lo eguaglio alle incognite alla destra

dell’uguale, lo faccio per più dati e quindi avrò tanti dati eguagliati con tante incognite. Ora vado a decidere

in che ordine mettere tutte queste equazioni, ad esempio in questo caso lo andiamo a fare procedendo per

colonne:

Prima prendo il primo elemento (in blu) poi tutti quelli sotto (in rosa), poi prendo il terzo elemento (in blu)

poi tutti quelli sotto (in giallo) e così via, un pezzo alla volta li ricavo tutti. Ora vado a concentrarmi sul primo

elemento:

i = 1 j = 1

Per la formula vista in precedenza in rosso avremo:

B B B B √

 

= = =

Ora ci andiamo a ricavare tutti gli elementi del diagrammino segnati come 2 in rosa:

i = 2,3,….,n j = 1

Quindi mi ricavo:

B B

= B

Ora e li conosco, quindi mi vado a ricavare:

B A

<

= AA 98

Ora ci facciamo un altro elemento:

Per un qualsiasi elemento la formula sarà:

Eccoci alla fine all’algoritmo di Cholesky:

Quindi abbiamo detto che Cholesky ci calcola la fattorizzazione e ci verifica che la matrice sia definita positiva.

La teoria ci dice che si verifica se il radicando, riportato qui sopra, viene > 0. Se strada facendo ci viene un

radicando nullo o negativo dovremo fermarci perché la matrice non è più definita positiva. Se prendiamo

Gauss va bene per matrici qualsiasi, ci dà la soluzione ma non è quella esatta quindi quando il problema è

simmetrico definito positivo, io devo usare Cholesky perché oltre a risolverli il problema correttamente mi fa

anche la verifica, ovvero se c’è un errore me lo dice. 99

Quindi se consideriamo di avere la mia matrice simmetrica e definita positiva e abbiamo:

T

A = L L A x = b

Prendo la fattorizzazione, la sostituisco al posto di A:

T

L L x = b

Poi vado a fare i soliti due passi:

• L y = b

• T

L x = y " !

D

La fattorizzazione mi costerà, essendo simmetrica, mentre il costo delle due triangolari sarà il solito di .

Quindi uso Cholesky, gli faccio fare la verifica, se c’è qualche problema mi viene segnalato, poi risolvo il

sistema. Matrici sparse

Questa è una piccola matrice presa da un problema ad elementi finiti, scannerizzata. E’ una matrice 500x500,

quindi piccolissima, in cui con i puntini sono evidenziati gli elementi non nulli, che sono 3971. Quando faccio

Cholesky, succede il cosiddetto fenomeno di fill-in (riempimento), questo succede anche con Gauss e con

tutti i metodi diretti. Quando abbiamo Gauss avremo:

/ E / E / E / E)

m(E

(E 0)

0 0 E 0 E 0 E

-

Quando vado a moltiplicare la prima riga per il moltiplicatore poi la vado a sottrarre a quella sotto, il risultato

verrà: (0 E)

0 E E E E E E

Quando vado a sottrarre, vado a mettere degli elementi non nulli dove prima cerano degli zeri, questo è il

fenomeno di fill-in, ed è intrinseco all’operazione che stiamo facendo. E’ un problema che per Gauss,

2

Cholesky e per i metodi diretti, devo tenere in memoria la mia matrice e predisporre n locazioni di memoria

perché non so come si riempiono. Ovviamente il riempimento dipende dalla disposizione dei non nulli. Quindi

in questo esempio vediamo che possiamo passare da 3971 elementi a 11533 elementi e questa cosa mi crea

un problema. 100

E’ possibile rimediare il problema fatto presente?

Si, andando a riordinare la matrice. Abbiamo dei metodi di riordinamento, in cui nelle immagini successive

ne vediamo l’esempio:

Vediamo che la matrice di prima con questo algoritmo assume questa forma, il totale è sempre lo stesso, è

solamente riordinata la matrice. Nel caso di Cholesky ora sono 9073 elementi. Ora vado a fare un altro

riordinamento:

La matrice è sempre la stessa mentre Cholesky ora ha 8440 elementi. Possiamo andare avanti e vedere quanti

altri me ne concede matlab e scegliere il nostro preferito.

Dei 3 quale sarà il nostro preferito?

Se siamo dei bravi ingegneri civili e facciamo dei buoni codici elementi finiti, la matrice deve essere il più

possibile a banda. Dei 3 il secondo sarà il nostro preferito perché la forma la so memorizzare bene e resta

tutta nella banda (ovvero ho tutti gli elementi non nulli nella diagonale). Quindi il secondo caso, anche se

apparentemente più costosa, però è la struttura che mi permette un’implementazione furba. Bisogna fare

attenzione perché i metodi di riordine sono costosi, quindi se la matrice viene usata più volte, allora posso

fare la fatica di fare un ordinamento, diversamente se la matrice è talmente grande che ho difficoltà con

Gauss e Cholesky, vado a usare i metodi iterativi. Da ricordarsi che la matrice è sparsa quindi ho pochi

elementi non nulli nelle nostre applicazioni. 101

Prodotto scalare

È una funzione che parte da due vettori di componenti reali o complesse e come risultato si ottiene un

reale o complesso. Tale funzione verifica le seguenti proprietà:

1. Se prendo un vettore qualsiasi e lo moltiplico scalarmente per sé stesso ottengo sempre un

risultato reale maggiore o uguale a 0. Si ottiene 0 se e solo se il vettore è quello nullo.

2. Proprietà di simmetria, in ambito reale x scalare y è uguale a y scalare x, in ambito complesso non

devo scordarmi del coniugato.

3. Proprietà di linearità, x,y,z, sono tre vettori qualsiasi, mentre α e β sono due scalari:

Esempi:

1 1

2 2 ⟨,⟩

=[ ] =[ ] = + + ⋯ +

1 1 2 2

⋮ ⋮

T

Lo stesso identico risultato lo si poteva ottenere moltiplicando il vettore y trasposto (y ) per x.

→

H

- se ho una matrice n x n scalare x Ax > 0 se A è definita positiva

→

H

- se ho una matrice n x n scalare x Ax ≥ 0 se A è semidefinita positiva

→

H

- se ho una matrice n x n scalare x Ax < 0 se A è definita negativa

→

H

- se ho una matrice n x n scalare x Ax ≤ 0 se A è semidefinita negativa

Con ciò vediamo che quando abbiamo una matrice che possiede la proprietà di simmetria ed è definita

positiva, possiamo definire un prodotto scalare.

nxn

Inoltre, se una matrice A è simmetrica definita positiva, anche tutte le sue sottomatrici principali sono

definite positive.

Teorema Criterio di Sylvester: serve per verificare se una matrice è definita positiva. Una matrice



simmetrica A di ordine n è definita positiva se e solo se det(A ) 0 k = 1,2,...,n ove A è la matrice principale

k k

di ordine k, cioè la matrice formata dalle prime k righe e colonne della matrice A.

102

Norme vettoriali

È una funzione che parte da un vettore di componenti reali o complesse, ma come risultato si ottiene

sempre un numero reale positivo più lo zero. Tale funzione verifica le seguenti proprietà:

1. Calcolare la norma di un vettore x qualsiasi il risultato è sempre un numero reale positivo, e viene 0

se e solo se x è il vettore nullo.

2. Se è presente uno scalare all’interno del vettore x, si può portare fuori con tanto di modulo.

3. Disuguaglianza triangolare, la norma della somma di due vettori è uguale alla somma delle norme.

Esistono infinite tipologie di norme, ma queste 3 sono le più utili:

- Norma 1: è la somma in valore assoluto delle sue componenti

- Norma infinito – Chebyshev: considera in valore assoluto tutte le componenti e poi prende la più

grande

- Norma 2 – euclidea: è la radice del prodotto scalare di x con sé stesso

Esempio:

La norma dal punto di vista numerico qual è la più conveniente? È la Norma 2 – euclidea perché è l’unica

che deriva da un prodotto scalare, perché ci definisce il concetto importante di ortogonalità; questo

concetto permette di semplificare gli algoritmi.

Esempio:

Si parte da dei dati sperimentali di un certo fenomeno. Si

deve costruire un modello continuo che descrive questo

fenomeno il più precisamente possibile.

Se il modello scosta di un certo valore dal punto che

rappresenta il dato sperimentale vuol dire che ci sarà un

certo errore. 103

Un buon modellino deve rendere minima la lunghezza del vettore errore.

Confrontiamo la norma 1 con la norma infinito e con la norma 2.

2 12 22 32 2

‖‖ = + + + ⋯ +

2

Se ad esempio il risultato fosse un numero piccolo (10) non si può dire che il modellino abbia tanti errori

piccoli, perché non posso sapere se sono uniformemente distribuiti gli errori, infatti posso avere

paradossalmente tutti errori pari a zero ed un unico errore pari a 10.

Per questo motivo la norma 2 non dà informazioni sulla distribuzione di errore e se la utilizzo accetto anche

il caso più estremo.

Se si considera la Norma infinito – Chebyshev, si ottiene come risultato il valore di norma della componente

massima, cioè la componente di picco di errore. Allora il nostro obiettivo è quello di minimizzare la norma

del picco massimo di errore perch&ea