Appunti completi del corso "Metodi analitici e numerici per l'ingegneria"

Errore di troncamento

Problema fisico/ingegneristico

Modello matematico

Approssimazione numerica

Tipi di errore

• Errore di arrotondamento: Soluzione esatta modello (roundoff) (lo trascuriamo)
• Soluzione approssimata
• Parole chiave: Approssimazione
• Errore assoluto:
• Errore relativo:
• Errore del modello: Errore che si ha nel passaggio tra il problema e il modello matematico, perché si introducono delle idealizzazioni.
• Errore di troncamento: Si sostituisce l'integrale (somma infinita), con una somma finita.
• Errore di arrotondamento

Software per soluzioni esatte

Esistono anche dei software che permettono di avere soluzioni esatte partendo dal modello matematico, come "Wolfram Alpha", quindi non si fanno approssimazioni.

Modello ben posto

La soluzione dipende con continuità dai dati, se si cambiano poco i dati, cambia poco la soluzione. (questo garantisce la validità delle soluzioni approssimate, che appunto si discostano un po' dalle soluzioni esatte)

Aritmetica esatta vs. Aritmetica della macchina

Roundoff

I computer non possono immagazzinare tutti i numeri, ma solo alcuni. Non è più vera la proprietà associativa: Floating Point normalizzati (con "normalizzati" si intende che la prima cifra dopo la virgola è ≠ 0)

Secondo questa trattazione però non si può scrivere lo zero. Questi però non sono tutti i numeri rappresentabili da un computer, ma sono solo quelli "normalizzati" rappresentabili.

Floating Point rappresentabili

Disegniamo sulla retta i numeri rappresentabili: i punti neri sono i floating points normalizzati. Fissato un dato esponente, vediamo cosa succede i punti sono tutti equidistanti se prendiamo:

I numeri rappresentabili sono pochissimi, e non sono equispaziati, sono molto fitti per numeri piccoli, e sono via via più radi per numeri più grandi (quindi non si perde tanto in errore relativo). Non si riescono però a rappresentare in numeri tra -0,25 e 0,25, per questo in questo intervallo si usano i:

Floating Point denormalizzati

I computer riempiono quell'intervallo con numeri (denormalizzati) con la stessa precisione degli intervalli adiacenti.

• L'errore assoluto sarà sempre lo stesso
• L'errore relativo sarà maggiore

Macchina epsilon (ε)

"Eps" restituisce la distanza di un punto dal successivo. Di fatto ε macchina è l'errore relativo: ε machina è importante perché permette al computer di approssimare un numero tra due punti, al punto più vicino, minimizzando quindi l'errore relativo. Quindi prendendo un numero qualsiasi, l'errore relativo che il computer commette nell'approssimazione è sempre inferiore a metà di ε machina.

Overflow e sistemi di equazioni lineari

Se x è maggiore di xmax non esiste la sua rappresentazione floating points, e siamo in un caso di overflow. Vediamo quanto valgono alcuni valori particolari:

MatLab:

• Proprietà che saltano:
• Proprietà associativa: (perché xmax + 1 fa ∞, meno xmax, che è un numero, fa sempre ∞)
• (perché gli si sta aggiungendo un numero troppo piccolo)
• (se x è molto piccolo, si sta dividendo per una cosa molto piccola)

Sistemi di equazioni lineari

Un sistema lineare contiene n polinomi di I grado con n incognite vettori colonna matrice Teorema esiste un'unica soluzione del sistema se, e solo se, la matrice A è NON singolare (det ≠ 0) Corollario la soluzione x è l'inversa di A moltiplicata per il vettore dei dati b.

Metodi "algebrici" per trovare:

• Calcolare l'inversa (al posto della i-esima colonna si mette il vettore b)
• Cramer (se n = 1)
• Metodo di Laplace per calcolare il determinante: Il metodo di Cramer è molto inefficiente, perché per calcolare il determinante sono necessarie molte operazioni elementari.

Metodo di eliminazione di Gauss

Per risolvere i sistemi lineari vedremo diversi algoritmi di soluzione:

• Sistemi lineari
• Metodi diretti (Analitici/Algebrici) (MEG o simili)
• Metodi iterativi (Numerici) (non danno una soluzione esatta in un numero finito di passi)

(Costo minimo possibile n operazioni, perché una matrice ha n componenti, e serve almeno un'operazione per coefficiente)

Metodi diretti: Metodo della fattorizzazione LU

Non singolare = matrice triangolare inferiore anche L e U devono essere non singolari = matrice triangolare superiore. A ammette una fattorizzazione LU, se A = LU. Supponiamo quindi che siamo passati da avere 1 sistema ad averne 2, ma sono molto più semplici.

Risoluzione di sistemi con matrici triangolari

• Il fatto che L e U non siano singolari, garantisce che il loro det sia ≠ 0, e siccome il det delle matrici triangolari si trova moltiplicando gli elementi sulla diagonale, ciò garantisce che l₁₁, l₂₂ e l₃₃ siano ≠ 0.

Algoritmo della sostituzione in avanti

Costo computazionale: (numero di operazioni necessarie) divisione, sottrazioni, moltiplicazioni.

Risoluzione di sistemi con matrici triangolari superiori

Algoritmo della sostituzione all'indietro. Il costo computazionale è uguale all'algoritmo precedente. Quindi il costo computazionale totale di questo metodo è:

Calcolo di L e U

Ma come possiamo calcolare L e U, ammesso che esistano? Dati incognite Uguagliare queste due matrici significa risolvere il sistema:

Abbiamo 4 equazioni e 6 incognite. Quindi avremo 2 GdL. In questi casi per convenzione gli elementi sulla diagonale di L vengono presi = 1.

• Equazioni incognite incognite
• Quindi ci sono n incognite in più, che però abbiamo visto che per convenzione si assumono uguali a 1, e sono proprio gli n elementi della diagonale della matrice L.

Algoritmo MEG (Eliminazione di Gauss)

(In realtà quello che vediamo adesso è il "metodo della fattorizzazione") Attraverso vari passaggi di sottrazione delle righe si possono costruire le matrici triangolari L e U. I moltiplicatori sono proprio gli elementi di L:

Introduzione dei moltiplicatori

Si introducono dei moltiplicatori: Gli elementi sulla diagonale sono anche chiamati "pivot", e devono essere ≠ 0. Si sottrae alla seconda riga la prima moltiplicata per il corrispettivo moltiplicatore: Si procede così anche per azzerare gli elementi della seconda colonna: moltiplicatore:

Quindi: L si ricava sapendo che sulla diagonale per convenzione ci sono tutti 1, e i suoi elementi sono proprio i moltiplicatori. Se ora si moltiplicano tra loro L e U, si ritrova proprio A.

Costo computazionale

Vogliamo adesso risolvere:

• Ipotesi:
• A sia singolare:
• Esista una fattorizzazione LU:

Ci sono tre modi per risolvere questo sistema:

• Gauss MEG (vero)
• Con l'inversa
• Matrice orlata Gauss

Costo: Quindi trovare è come risolvere n sistemi. Conviene usare il primo metodo anche se ha un costo maggiore, perché se per esempio si ha un caso del genere: Gauss fatto su A basta farlo una volta, e vale per tutte le q, invece con il secondo metodo bisogna applicare di nuovo Gauss ad ogni passaggio.

Condizioni per applicare il metodo di Gauss

Non è sempre possibile applicare il metodo di Gauss: (bisognerebbe prima verificare che det(A) ≠ 0) l'algoritmo di Gauss non funziona. Come si poteva capire prima? (condizione necessaria e/o sufficiente per di una fattorizzazione) Se det(A) ≠ 0 ma A ≠ LU cosa si può fare?

Per la prima domanda troveremo:

• 1 condizione N&S
• 3 condizioni S

Condizione N&S

Teorema: se le sottomatrici principali di hanno tutte det ≠ 0 allora regola Nord-Ovest:

Condizioni S

• A matrice a "dominanza diagonale stretta per righe" (preso un qualunque elemento sulla diagonale, questo elemento domina la riga)
• A matrice a "dominanza diagonale stretta per colonne"
• A simmetrica, definita positiva (sdp) prodotto scalare autovalori

A è simmetrica, definita positiva, se e solo se, tutti i suoi autovalori sono positivi: (Se una matrice è simmetrica, tutti i suoi autovalori sono reali. Una matrice simmetrica è definita positiva se e solo se tutti i suoi autovalori sono positivi) Se una di queste è verificata, si può rispondere alla prima domanda. Per matrici quadrate vale che: autovettore autovalore Pivoting Risposta alla seconda domanda: tecnica di "pivoting" (scambio di righe)

Pivoting

Se si arriva ad un certo punto (passo i-esimo) dove un pivot è zero, si scambia la i-esima riga con la successiva. Lo scambio di due righe si chiama: "pivoting per righe". Per tenere traccia di questo scambio si usano le "matrici di permutazione". Sono matrici ottenute dalla matrice identità scambiando righe (o colonne). ! bisogna moltiplicarla a sinistra !

Oss: se il pivot ma è molto piccolo, l'algoritmo prosegue ma con grossi errori. Non bisogna per forza fare lo scambio con la riga immediatamente successiva, ma si può fare anche con le alte righe successive. Di solito conviene scegliere di scambiare la riga con il pivot più grande.

Se scambio riga i con riga k > i, dove in realtà conviene sempre applicare questo ragionamento, cioè fare sempre lo scambio, anche se i pivot sono ≠ 0, ma sono piccoli ( = non sono il più grande). Quindi non si fa lo scambio solo se il pivot a quel passaggio è il più grande possibile.

Quindi il pivoting viene fatto per:

• Per necessità (aii = 0)
• Per accuratezza (aii ≠ 0 , ma piccolo)

"Pivoting per righe": ad ogni step scambio 2 righe per massimizzare il modulo del pivot. Gli scambi tra righe effettuati su A, vengono replicati anche su P. Questo perché quando si arriva in fondo si trova una fattorizzazione LU, ma non è una fattorizzazione della matrice A, ma della matrice A a cui sono state scambiate delle righe, che è uguale proprio a PA. Alla fine:

Accuratezza del pivoting

Vediamo perché bisogna usare il pivoting anche per "accuratezza": A è fattorizzabile senza pivoting LU (2) Il problema è che si avrà un pivot a molto piccolo, il che genera dei moltiplicatori molto grossi. Se si facesse A - LU in teoria si dovrebbe avere una matrice con tutti zeri: è dovuto alla propagazione dell'errore di questo piccolo decimale. Quindi la fattorizzazione NON è accurata. Invece la fattorizzazione con pivoting è accurata, perché permette di scambiare le righe in modo da non avere il pivot piccolo che generava inaccuratezza. Essere accurata significa che PA - LU (che dovrebbe avere tutti zeri), ha comunque termini molto piccoli.

Comando MatLab per la fattorizzazione: [L,U,P] = lu(A)

Algoritmi per particolari A

Esistono algoritmi più economici per particolari A:

• Se A simmetrica, definita positiva: "Algoritmo di Cholesky" Si trova una fattorizzazione di questo genere: estensione di questa. Dato che le matrici L e U sono una la trasposta dell'altra, bisogna fare la metà dei calcoli. Proprietà (in questo caso non si prendono gli 1 sulla diagonale di L, ma sulle diagonali ci sono le radici degli autovalori positivi)
• Se A è tridiagonale: "Algoritmo di Thomas" Questo algoritmo è rigorosamente senza pivoting, perché se si scambiassero delle righe o delle colonne, la matrice non sarebbe più tridiagonale.

Sostanzialmente è sempre l'algoritmo di Gauss, ma il fatto che ci siano tanti zeri semplifica molto i calcoli. Il costo computazionale di questo algoritmo infatti è:

Tutti i problemi in cui si prende una struttura continua e la si discretizza, tipicamente sono problemi di questo genere, dove le matrici corrispondenti sono:

• "Matrici sparse": il numero di termini diversi da 0 è dell'ordine di n: # { termini ≠ 0 } = o(n)
• "Matrici strutturate": i termini diversi da 0 sono posizionati bene: { termini ≠ 0 } "posizionati bene"

Errori e accuratezza

Quanto è accurata la soluzione di ? "Analisi di sensibilità (sensibilità)"

Fattori di accuratezza

Cose che abbiamo fatto:

1. Fattorizzazione LU con pivoting: se è accurata vuol dire che (e anche l'errore relativo)
2. Calcolato la soluzione (la soluzione vera però è ): se è accurata vuol dire che

Chiediamoci adesso: se la fattorizzazione LU è accurata, il calcolo della soluzione è accurato? In generale NO Prendiamo una matrice particolare, detta "Matrice di Hilbert": Consideriamo il sistema:

• Scegliamo in modo tale che (la soluzione quindi sarà un vettore con tutti 1)
• Calcoliamo la soluzione con MatLab: (questa soluzione è diversa da )

Vediamo qual è l'errore relativo: Vediamo anche l'errore nella fattorizzazione: è una matrice, consideriamo il maggiore tra i suoi elementi: Quindi anche se la fattorizzazione LU è accurata, la soluzione può essere anche molto sbagliata.

Norma

Per misurare gli errori, introduciamo il concetto di "norma".

Norme di vettori

• Proprietà di base:
• Definita la norma, si possono definire anche i concetti di:
• "Distanza" tra due vettori:
• "Convergenza": Se non si specifica, ci si riferisce alla "norma euclidea"
• Palla unitaria (punti a distanza 1 dall'origine):

Altre norme

Punti raggiungibili percorrendo una distanza pari a 1, e muovendosi solo lungo x e y (utile per calcolare il percorso in una città, in cui si devono seguire le strade, non si può tagliare per le case).

Sono tutte equivalenti, tutte convergono allo stesso modo. Sono solo modi diversi per misurare le distanze.

Norme di matrici (quadrate)

Noi vediamo una famiglia di norme dette "norme operatoriali", che soddisfano sempre le proprietà e prima.

Norma spettrale

È la controparte della norma euclidea: è fatta in modo tale che: (la norma del vettore Ax è <= a una costante per la norma di x) (è un vettore) è la più piccola possibile norma euclidea la matrice A raddoppia le coordinate orizzontali, e lascia invariate quelle verticali. Quanto può aumentare o diminuire la lunghezza di x?

Sicuramente l'aumento è più piccolo di 100 volte, ma anche di 10 volte, ma qual è la costante più piccola?

Se la matrice è simmetrica, ammette base ortonormale di autovettori, cioè è diagonalizzabile (teorema spettrale) autovalore. Il metodo è lo stesso del caso precedente, usando però come assi gli autovettori. Autovettore se non è definita positiva se A è simmetrica, definita positiva la norma è il più grosso degli autovalori.

Definizione generale

Definizione: sia qualunque. Se è simmetrica quindi torniamo al caso di prima -1. Se A è invertibile (det ≠ 0), gli autovalori di A sono 1 / gli autovalori di A: dato che c'è l'1 fratto, dobbiamo prendere il più piccolo autovalore o è perché gli autovalori dell'inversa sono 1/..., e qua si stanno prendendo gli autovalori di A normale???

Risolvere numericamente

Esatta (reale sistema risolto dalla macchina) macchina Errore: Errore relativo: Residuò: (idealmente dovrebbe essere 0, se la soluzione fosse quella esatta) Residuò relativo: ! senza conoscere non si possono calcolare gli errori, ma si possono calcolare i residui ! ma se r è "piccolo", e è "piccolo"?

Rel e rel Caso più semplice (ma verosimile usando pivoting): è l'errore è il residuo. Introduciamo le norme: Moltiplichiamo (1) e (2), si possono moltiplicare membro a membro perché i numeri sono tutti positivi: numero di condizionamento (spettrale) della matrice. Passando dal residuo relativo, all'errore relativo, si perdono tanti ordini di grandezza, quanti sono quelli del numero di condizionamento della matrice.

"Matrice ben condizionata" "Matrice mal condizionata" Esempio: matrice di Hilbert con n = 4:

Osservazioni

Oss. 0: se A è simmetrica, definita positiva se autovalore massimo e autovalore minimo sono uguali, A è un multiplo della matrice identità.

Oss. 1: se Oss. 2: MatLab:

- cond(A) = cond(A,2) Se non è specificato un diverso numero, restituisce il numero di condizionamento spettrale di A. Questo comando però non vale per matrici sparse.

- condest(A) Restituisce solo una stima di K (A), cioè il prodotto delle norme 1, non delle norme 2.

Norma 2 vs Raggio spettrale

Norma 2

La norma 2 di una matrice è la più piccola costante tale che: Affinché la norma 2 sia la più piccola possibile, bisogna prendere il più grande valore del rapporto ||Ax||/||x||:

Raggio spettrale

Il raggio spettrale è il massimo modulo degli autovalori. In generale norma spettrale e raggio spettrale non coincidono. Di sicuro si può dire che il raggio spettrale sarà sempre più piccolo della norma 2:

Possiamo affermare ciò perché se prendiamo un autovalore e un autovettore associato all'autovalore: Da questo confronto si vede che la norma 2 è maggiore o uguale del modulo di qualunque autovalore.

Proposizione 1

Se è una matrice sdp: (senza modulo perché tanto i λ sono positivi) Quindi se è u