Calcolo Numerico

Introduzione

Dato un processo fisico, è spesso possibile introdurre un opportuno problema matematico per la sua descrizione quantitativa. Un processo fisico può essere per esempio il sangue che scorre in un'arteria: questo fenomeno può essere descritto da un problema matematico. La scelta di questo problema è unica o ci sono tante scelte? La scelta del modello matematico non è unica, ma dipende dalle ipotesi di lavoro che faccio. Prima di scrivere il problema matematico è necessario fare ipotesi: è impossibile descrivere il sangue esattamente com'è, con globuli rossi, piastrine, rugosità della parete… più le ipotesi sono vicine alla realtà, più la descrizione del fenomeno con il modello matematico sarà accurata. Il prezzo da pagare sarà la complessità del modello. Ipotesi più rigide semplificano il modello perdendo in accuratezza. Ipotesi modellistiche lasche -> Grande accuratezza e complessità La scelta del modello dipende dalla capacità di risolvere un modello complesso, il tempo a disposizione per risolverlo… è la quantità di interesse fisica, quella cosa che vorrei conoscere ma non conoscerò mai.

#d sono i dati (per esempio le caratteristiche del singolo paziente) X è la soluzione continua matematica del problema. In generale non è possibile risolvere analiticamente i problemi. È necessario introdurre metodi numerici per la soluzione approssimata di X = 0. I metodi numerici permettono di risolvere in maniera approssimata il problema continuo. h è la soluzione numerica del problema.

In che relazione sta X con h, e in che relazione sta h con X?

Ho quindi due fonti di errore:

• Errore di modello: è l'errore commesso pretendendo di descrivere la realtà con la matematica. X = X - h
• Errore numerico: X = X - h

Idealmente il metodo numerico deve soddisfare due requisiti:

• Rendere l'errore numerico piccolo
• Produrre una soluzione h in un tempo computazionale accettabile

Così come ci sono tanti problemi continui descriventi lo stesso processo fisico, dato un problema continuo è possibile introdurre diversi problemi numerici.

Risoluzione di sistemi lineari

3x3 A ∈ R, b ∈ R^h < A, b ∋ In forma per componenti: 3; A(x) = b, x = 1, ..., <= A = <=>1 In realtà abbiamo una formula analitica per risolvere un sistema lineare: è la formula di Kramer: det A < = det BCCD Dal punto di vista del costo computazionale, un metodo numerico consiste di fatto in una sequenza di operazioni elementari (+, −, ×, :): questo è un algoritmo. Questo va implementato al calcolatore, che svolge velocemente le operazioni elementari. È importante quantificare a priori, se possibile, il tempo di calcolo Anche per la formula di Kramer, nonostante sia una formula analitica, esiste un algoritmo: è possibile passarlo direttamente al calcolatore come algoritmo. Qual è il costo computazionale della formula analitica?

Se l'operazione C costa asintoticamente (ossia per m che tende a infinito) m operazioni elementari, scrivo ~. ~( )! ~( )( )! det A + 1 → + 1 + 11

Ipotizzo un calcolatore che esegua 10 flops/s = 15¹³ 14( )!+ 1 ≈ 2 ⋅ 10 → ≈ 10 → 20 Immagino che il sistema diventi un po' più complesso: A = 20¹⁹ 20 9( )!+ 1 ≈ 10 → ≈ 10 → 10 A = 100 A = 100 160( )!+ 1 ≈ 10.

Casi particolari

Considero dei casi particolari, in cui il sistema lineare è risolvibile bypassando la formula di Kramer.

Matrice diagonale

È una matrice che ha zeri ovunque tranne che sulla diagonale:

\ 0U W10 \ = 0 ≠ Se la matrice è diagonale, le equazioni del sistema sono disaccoppiate tra loro. A = 3 <; b = A → x = A → A =<= A = < << < < < <==>1 Questo sistema è risolvibile con un singolo ciclo for: il costo complessivo è n.

Matrice triangolare inferiore

A = 0 < \ 0U W \3 1 = 1: b = A → x = A → A =1= A = 1 11 1 1 1 11=>13 − 2 21 1 = 2: b = A → A + A → A =2= A = 2 21 1 22 2 2 22=>1 ∑ − A < =]< <= A = A << Dal punto di vista di codice, c'è un ciclo for che scorre la i, e un secondo ciclo for annidato che esegua la sommatoria su j. Ci aspettiamo che il costo computazionale sia asintoticamente n². Questo metodo è detto metodo delle sostituzioni in avanti.

Matrice triangolare superiore

È analogo al caso della matrice triangolare inferiore.

Caso generale

Def: A, i = 1, ..., n sono le sottomatrici principali di A e sono ottenute prendendo le prime i righe e le prime i colonne di A. Teo: se i (A) ≠ 0, i = 1, ..., n - 1, allora esistono un'unica matrice triangolare inferiore L con L = 1 ed un'unica matrice triangolare superiore U, tali che A = LU. Questo metodo è detto fattorizzazione LU. Questa fattorizzazione è unica, a patto che gli elementi della matrice L siano fissati uguali a 1. Supponiamo che esista la fattorizzazione LU. A² = ~A = LU ↔ A = LU ↔ a →A² = ~A². Il costo A è il costo della risoluzione dei sistemi. Ma quanto mi è costato trovare L e U? Devo capire come e quanto costa trovare L e U. Per trovare i fattori L e U si usa l'algoritmo dell'eliminazione gaussiana (MEG), che è basato su n passi, ad ognuno dei quali aggiorno tutti gli elementi della matrice.

Il costo sarà A³, perché ciascun passo ha bisogno di 2 cicli for annidati per aggiornare ogni singolo elemento della matrice. Il metodo della fattorizzazione LU nella sua interità consta A³. Una matrice 100.000X100.000 costa 10 secondi, ossia circa 3 ore. Fino ad una certa dimensione delle matrici è piuttosto accettabile, oltre un certo livello è comunque proibitivo. L'idea alla base dell'algoritmo è che ad ogni passo k, si elimina l'incognita A dalle righe k + 1, ..., n. Alla fine, la matrice sarà la triangolare superiore U. Il passaggio A = A è nascosto negli n − 1 passi da k = 1, ..., n − 1. l'ultimo passo è k = n.

Esempio

1 1 1 3l m = 2 3 57 8 93. Elimino A dalla II e dalla III riga. Sostituisco in modo da avere 0 nei primi due termini della seconda e della terza riga. 1 1 3 → b − 2 pn l m→ = 0 1 −1 → b − 7 0 1 −12p4. Elimino A dalla III riga di 2 1 1 3p l m → b − b→ A = = 0 1 −10 0 −115. Risolvo A = La y è il frutto delle operazioni somma-sottrazione fatte su b. per esempio5 1q r = −23= 1 1 = A − 2 = −4 t2 2 1s q r→ A = −4 = −7 = −43 3 2 −4= 21q r = −40 Supponiamo di avere un sistema lineare con la stessa matrice A, ma con termine noto b diverso. Cosa cambia nella procedura? Cambia solo y, quindi se è vero il teorema di fattorizzazione unica la U non cambia. La fattorizzazione di una matrice A può essere quindi fatta una volta sola, a prescindere dalla matrice b del problema. Notiamo che data A, la fattorizzazione LU è unica, quindi il costo A è fatto una tantum: se cambio A non devo ricalcolare L e U.

Esempio 2

3 6 9l m = 2 4 57 8 91. Sostituisco: 2 → b − 3 6 93 p l mv → = 0 0 −17 0 6 −12 → b − 3 Questo è un problema! Non voglio zeri sulla diagonale principale 3 6 9l m2. ... A = 0 0 −1− − − 3. Non posso fare il terzo passo! Data una matrice triangolare, gli autovalori sono A = b → b , ∃ b = 0 → det A = 0 → A = A è << < e il metodo di eliminazione gaussiana si interrompe. È un difetto a monte, ossia il problema continuo non è risolvibile, o il problema è risolvibile ma la tecnica numerica in questo caso fallisce? Per rispondere a questa domanda devo calcolare il determinante di A. det A = 18 ≠ 0 In questo caso il sistema lineare A = b è risolvibile, ma il metodo dell'eliminazione gaussiana fallisce. Perché? In questo esempio non è soddisfatta l'ulteriore ipotesi necessaria per applicare LU: è necessario che anche i determinanti delle sottomatrici principali siano non nulli, che in questo caso non è vero. det A = 12 − 12 = 02 p In generale, se det A = 0, allora al passo k = 0 → A = 0 e il sistema non è risolvibile con MEG. zz

È necessario trovare una soluzione al problema dell'esempio 2! È sufficiente scambiare la II e la terza riga di A: a patto che si faccia la stessa cosa su b, la soluzione sarà identica. 6 3 6 9l m = 7 8 92 4 57 3 6 9 → b − p3{ l m4. → A = 0 −6 −122 → b − A 0 0 −13 p Questa tecnica è detta pivoting: se det A = 0, scambio righe di A al passo k per ovviare al fatto che p det A = 0. In questo modo evito A = 0. Non posso scambiare la riga k con una riga sopra di essa, perché al passo k le righe superiori sono già state modificate dall'algoritmo e non vanno più modificate. p Quindi scambio la riga k con una delle righe sotto (>k). Scambio con la riga l>k, purché A ≠ 0. Può essere necessario eseguire più volte il pivoting! Algebricamente il pivoting coincide con il moltiplicare la matrice A con una matrice P: A = PA. A = PA → PA = L → L = U ↔ | L = U

Questo metodo è numerico? Il MEG non effettua approssimazioni numeriche, quindi non sarebbe un metodo numerico. Ma lo è! In aritmetica esatta, il MEG non introduce alcun errore numerico. L'aritmetica esatta è il "mondo" in cui le operazioni elementari sono svolte esattamente. Questo mondo è la mente umana, che è priva di errori intrinsechi ed è in grado di svolgere esattamente le operazioni matematiche elementari. Questo non succede però al calcolatore, non per un suo limite operativo (non sbaglia a svolgere i calcoli!) ma per problemi di memorizzazione. Il calcolatore si trova quindi nel "mondo" dell'aritmetica finita. Questo introduce un errore nell'algoritmo. In generale un numero (come π, non periodico e con un numero di cifre decimali infinite) avrebbe bisogno di infinite celle di memoria. Se devo sommare due numeri π e e, di fatto sommo € e π, cioè la loro memorizzazione su un numero finito di celle: € = € + e. Questo errore di per sé è molto piccolo, svolgendo un numero elevato di operazioni si può accumulare un errore grande (errore di arrotondamento) che si propaga: ogni volta che si svolge una delle operazioni si introduce un piccolo errore che si accumula ai precedenti. In aritmetica esatta il MEG è un metodo esatto, ma di fatto quando è dato in pasto al calcolatore diventa un metodo numerico perché quello che il calcolatore restituisce è affetto da errore. Il MEG non ha errore numerico, ma ha un errore di algoritmo! la matrice che viene memorizzata dal calcolatore è la matrice effettiva sommata ad un piccolo errore.

• → A = A + | „ → A = A + (A)

Di fatto sto risolvendo A + A = A + ††. La speranza è che A = A + A, con A piccolo. La norma di una matrice è un numero scalare associato alla matrice che ne quantifica le dimensioni.

• • • „‡ La fattorizzazione LU è fatta su → A = LU.

Le operazioni dell'algoritmo sono svolte esattamente ma su numeri sbagliati. &per mil

Condizionamento

Teo. Sia A = LU. Se A è simmetrica, definita positiva (ossia A ⋅ x ⋅ x > 0 ∀ x > 0) allora A = LU ed è detto numero di condizionamento. Questo numero stabilisce quanto gli autovalori sono lontani tra loro: se k è elevato gli autovalori sono dislocati in un'ampia regione, se k invece è piccola gli autovalori sono confinati in una regione piccola. Ogni matrice ha il suo numero di condizionamento, che non può cambiare.

†’“† †’•†(A) Supponendo A = 0 (per esempio A è composta solo da numeri interi), si ha che ≤ A. †“† †•† Matrici intrinsecamente malcondizionate (con autovalori sparsi un po' ovunque) producono un errore di arrotondamento che si propaga molto: la soluzione effettiva A potrebbe essere molto lontana dalla soluzione corretta A. È quindi necessario stimare il numero di condizionamento, per sapere se è possibile applicare il MEG senza introdurre un errore troppo grande. Supponiamo di essere in una situazione con A grande. Dati A e A (€ e A ̂) la più sensibile fra le operazioni elementari è A - A. In generale vorrei A il più piccolo possibile: più A è grande, più l'errore di propagazione si propaga. Nel MEG al passo k divido per A e sottraggo questo numero ad un altro: A A; Proprio per questa divisione se A = 0 l'algoritmo si interrompe! Allora vorrei A ≠ 0.

—t ˜˜ p Vorrei inoltre che sia grande, in modo da propagare l'errore il meno possibile. Durante il pivoting posso scegliere di scambiare la riga k con la riga l>k, purché A ≠ 0. p p In particolare, prendo l tale che A = max , in modo da ridurre al minimo l'errore di arrotondamento. È quindi conveniente eseguire il pivoting ad ogni passo, anche se A ≠ 0. In questo modo la propagazione dell'errore si abbassa. Questo produce A < A.

Memorizzazione delle matrici

Un aspetto importante a livello pratico è la memorizzazione delle matrici. „ ‡ A = 1: non è necessario memorizzare tutti questi elementi. Di conseguenza per memorizzare e A ho bisogno esattamente di n² celle (la parte triangolare di U e la sua diagonale, e la parte triangolare inferiore di L senza la diagonale):

U W • •„ ‡ Non è necessario tenere in memoria A. Memorizzo e A sullo spazio che avevo dedicato ad A. Considero una matrice A sparsa, cioè con un numero di elementi nulli asintotico ad n (il numero degli elementi non nulli è asintotico ad n). Esempio \ \ 0\ \ l m: Matrice triangolare gli elementi non nulli sono solo quelli sulle tre diagonali principali. 0 \ \# A = A + n - 1 + n - 1 ∼ A. Raddoppiando n il numero di elementi della matrice quadruplica, mentre il numero di elementi non nulli raddoppia seguendo n, non n². n² # = n − 3 + 2 ∼ Le matrici sparse sono interessanti perché sono molto ricorrenti risolvendo equazioni differenziali, metodi degli elementi finiti… gli elementi memorizzati sono solo quelli non nulli, e la loro ubicazione. Di fatto l'occupazione di memoria è ~n celle. La matrice è memorizzata con il costo in memoria di un vettore. Il problema è il cosiddetto fill-in: se A è sparsa (non metto i cappucci, quindi stiamo considerando l'aritmetica esatta), allora in generale L e U non lo sono. Non è quindi possibile sfruttare l'occupazione di memoria utilizzata per A per LU, ma è necessaria più memoria. Per questa ragione la fattorizzazione LU è poco utilizzata per matrici sparse. Vorrei trovare un metodo che sfrutti la sparsità della matrice e non costringa ad utilizzare una memorizzazione matriciale. Non è vero che per il fill-in il risultato sarà sbagliato: è tutta una questione di memoria! Risolvere il sistema con un numero di condizionamento alto è vietato: il risultato potrebbe essere molto lontano dalla...