Appunti Metodi numerici

1 Slides 1

Avvalendosi di un calcolatore si lavora con una aritmetica finita che dà luogo ad errori.

• Errore inerente: Si commettono errori rappresentando i dati del problema con un numero finito di cifre.
• Errore algoritmico: Si commettono errori eseguendo le operazioni aritmetiche richieste dal processo di calcolo con una aritmetica finita.

Un problema è malcondizionato se piccole variazioni sui dati inducono grandi variazioni sui risultati.

Definizione: la complessità di calcolo, detta anche complessità computazionale o costo computazionale, viene misurata in funzione delle dimensioni del problema che si sta trattando e dà indicazione del numero di operazioni aritmetiche che necessitano per ottenere il risultato.

2 Slides 2

2.1 Rappresentazione floating point

q·x = p N.

Dove: p ∈ R è un numero reale; N è la base del sistema di numerazione scelto; q è un intero positivo → p mantissa, q esponente.

Esempio: x = 0.12345 = 12345 · 10⁻⁵.

Dove “12345” è la mantissa (p) e “-5” è l’esponente (q). Quindi ad x si associa la coppia (p,q).

La rappresentazione normalizzata, cioè −1 ≤ |p|N < 1, conduce alla rappresentazione unica di un numero x = 0 mediante la coppia (p, q).

Lo spazio riservato alla rappresentazione di un numero reale può essere visto come: segno (s), esponente (q), mantissa (|p|).

Si avrà che: per la singola precisione ci saranno a disposizione 32 bit quindi: 1 bit per il segno, 8 bit per l’esponente e 23 bit per la mantissa.

Per la doppia precisione ci saranno a disposizione 64 bit quindi: 1 bit per il segno, 11 bit per l’esponente e 52 bit per la mantissa.

−127 ≤ q ≤ 128 (semplice precisione).

−1023 ≤ q ≤ 1024 (doppia precisione).

Definizione: x ∈ R è detto numero macchina.

q ∈ / m, M. Se q < m si associa zero ad x ed il sistema segnala la situazione di underflow.

Se q > m, x non viene rappresentato ed il sistema segna la situazione di overflow.

q q ∈ ±pN ∈ ℑ.

Se q ∈ [m, M] ma x = / allora si usa l’approssimazione x = p N_e e_e.

2.2 Troncamento o chopping

Nella mantissa p vengono escluse le cifre dopo la t-esima.

p = 0.a₁ a₂ a₃ ... a_t a_t+1...

p_e = 0.a₁ a₂ a₃ ... a_t.

Errore assoluto: ε_A = x_e − x.

Errore relativo: ε_R = (x_e − x) / x.

|ε_A| < N^−t+q.

|ε_R| < N^−t+1.

2.3 Arrotondamento

x_e = arr(x), x_e = p_eN^q.

trn(p) se a_t+1 < N/2.

trn(p + N^−t) se a_t+1 ≥ N/2.

|ε_A| < 1/2 N^−t+q.

|ε_R| < 1/2 N^−t+1.

Il valore che limita superiormente l’errore relativo (eps) è detto precisione di macchina.

eps = N^−t+1 troncamento.

eps = 1/2 N^−t+1 arrotondamento.

2.4 Operazioni di macchina

Se si effettua un’operazione su due numeri di macchina, il risultato non per forza lo è. Bisogna approssimare affinché: x, y ∈ ℑ → x op y / ∈ ℑ ma fl(x op y) ∈ ℑ.

ε_R = fl(x op y) − x op y / x op y.

|ε_R| < eps.

Definizione: fl(x op y) = (1 + ε_R)(x op y).

In aritmetica non è detto che valgano le proprietà dell’aritmetica:

fl(fl(x · y) · z) ≠ fl(x · fl(y · z)).

fl(fl(x + y) + z) ≠ fl(x + fl(y + z)).

Esempi da pagina 18 Slides 2.

2.5 Errore inerente

Errore relativo alla rappresentazione.

ε_In = f(x_e) − f(x) / f(x).

Sviluppo in serie di Taylor al primo ordine f(x):

f(x_e) = f(x) + δf(x)/δx_i (x_ei − x_i) + ... + δf(x)/δx_n (x_en − x_n).

f(x_e) − f(x) = δf(x)/δx_i (x_ei − x_i) + ... + δf(x)/δx_n (x_en − x_n).

Divido entrambi i membri per f(x) e moltiplico e divido il secondo membro per x_i ottenendo:

ε_IN = Σ_i=1ⁿ (x_ei − x_i) / x_i · x_i / f(x) · δf(x)/δx_i.

Sapendo che ε_ri = (x_ei − x_i) / x_i si può scrivere:

ε_IN = Σ_i=1ⁿ ε_ri x_i / f(x) · δf(x)/δx_i.

Ponendo C_i = x_i / f(x) · δf(x)/δx_i si ottiene:

ε_IN = Σ_i=1ⁿ ε_ri C_i.

Dove C_i è l’indice di condizionamento.

2.6 Errore algoritmico

ε_ALG = Σ_i=1^M θ_i (x₁, x₂, x₃, ..., x_n) γ_i.

Dove M è il numero di operazioni che compongono l’algoritmo in uso.

γ_i := ε_Ri, i = 1, 2, 3, ..., M.

2.7 Errore totale

ε_TOT ≈ ε_IN + ε_ALG.

Esempio da pagina 33 Slides 2.

3 Slides 3

a₁₁x₁ + a₁₂x₂ + ... + a_1nx_n = b₁.

a₂₁x₁ + a₂₂x₂ + ... + a_2nx_n = b₂.

........................................................................................

a_m1x₁ + a_m2x₂ + ... + a_mnx_n = b_m.

Da cui si può ottenere la matrice A associata al sistema, il vettore x delle incognite, e il vettore b dei termini noti. Quindi si può scrivere il sistema come: Ax = b.

Esempi da pagina 4 Slides 3.

3.1 Norma di un vettore

La norma di un vettore è un’applicazione che a un vettore associa un numero reale.

|| · || : Cⁿ −→ R⁺₀.

Definizione:

• 1. ||x|| ≥ 0 ∀x ∈ Cⁿ, ||x|| = 0 ↔ x = 0.
• 2. ||αx|| = |α| ||x|| ∀x ∈ Cⁿ, ∀α ∈ C.
• 3. ||x + y|| ≤ ||x|| + ||y|| ∀x, y ∈ Cⁿ.

Le funzioni norma più comunemente usate sono:

||x||₁ = Σ_i=1ⁿ |x_i|, norma 1.

||x||_∞ = max_{1 ≤ i ≤ n} |x_i|, norma infinito.

||x||₂ = √(Σ_i=1ⁿ |x_i|²), norma 2.

Algoritmo pagina 12 Slides 3.

3.2 Norma di una matrice

La norma di una matrice è un’applicazione che a una matrice associa un numero reale.

|| · || : C^n×n −→ R⁺₀.

Definizione:

• 1. ||A|| ≥ 0 ∀A ∈ C^n×n, ||A|| = 0 ↔ x = Ω.
• 2. ||αA|| = |α| ||A|| ∀A ∈ C^n×n, ∀α ∈ C.
• 3. ||A + B|| ≤ ||A|| + ||B|| ∀A, B ∈ C^n×n.

Le funzioni norma più comunemente usate sono:

||A||₁ = max_{1 ≤ j ≤ n} Σ_i=1ⁿ |a_ij|, norma 1.

||A||_∞ = max_{1 ≤ i ≤ n} Σ_j=1ⁿ |a_ij|, norma infinito.

||A||₂ = √(ρ(A^TA)), norma 2.

Dove ρ è il raggio spettrale cioè il massimo tra gli autovalori in modulo.

Algoritmo pagina 12 Slides 3.

3.3 Relazioni tra norma di vettore e di matrice

Compatibilità tra norma di vettore e di matrice:

||Ax||_∞ ≤ ||A||_∞ ||x||_∞.

Norma di vettore 1, 2 e ∞ sono compatibili con una norma 1, 2 e ∞ di una matrice.

3.4 Condizionamento di un sistema lineare

K(A) = ||A|| ||A⁻¹||.

In MATLAB, cond(A) fornisce in output l’indice di condizionamento del sistema Ax = b.

Da pagina 19 Slides 3.

3.5 Metodi numerici per la risoluzione di sistemi di equazioni lineari

Da pagina 22 Slides 3.

3.6 Metodi diretti per la risoluzione di sistemi lineari

Sono dei metodi basati sulla trasformazione del sistema in uno equivalente che ha una struttura più semplice, e quindi è più facile calcolare la soluzione in un numero finito di passi.

Vengono usati per sistemi di dimensioni non troppo grandi e con la matrice del sistema piena, cioè con la maggior parte degli elementi diversi da zero.

3.7 Metodo di eliminazione di Gauss (MEG)

Il Metodo di eliminazione di Gauss trasforma in n-1 passi il sistema lineare A x = b, A ∈ R^n×n, x, b ∈ R^n×1, det(A) ≠ 0, in una matrice dei coefficienti triangolare superiore. Ha un costo computazionale O(n³/3).

Formule generali:

Per k = 1, 2, ..., n-1 e i = k+1, k+2, ..., n.

m_ik = a_ik^(k) / a_kk^(k).

a_ij^(k+1) = a_ij^(k) − m_ik a_kj^(k).

b_i^(k+1) = b_i^(k) − m_ik b_k^(k).

Esempio da pagine 16 Slides 4.

Algoritmo Gauss (prima versione).

Input: A, n.

Output: Aⁿ, bⁿ.

• Ripetere per k = 1, n-1.
• Ripetere per i = k+1, n.
• m_ik = a_ik / a_kk.
• Ripetere per j = k+1, n.
• a_ij = a_ij − m_ik a_kj.
• b_i = b_i − m_ik b_k.

3.8 Metodo di Gauss con pivoting parziale

Se a_kk^(k) = 0 si cerca un r ≥ k tale che |a_rk^(k)| = max_k≤s≤n |a_sk^(k)|.

Cioè si scambia la riga k-esima con la riga r-esima. Conviene, ad ogni passo k, usare questa tecnica, anche se a_kk^(k) ≠ 0 perché aumenta la stabilità numerica.

Esempio da pagina 34 Slides 4. Codice pagina 36 Slides 4.

4 Slides 4

4.1 Fattorizzazioni

Il sistema Ax = B viene risolto operando una fattorizzazione della matrice A = BC.

Questo significa che BCx = b e cioè che: Cx = y e By = b.

4.2 Fattorizzazione LU

La scelta delle matrici B e C caratterizza la fattorizzazione LU. Avremo che:

B = L, L è triangolare inferiore.

C = U, U è triangolare superiore.

A = LU.

La matrice A, con elementi a₁₁, a₁₂, a₁₃, ..., a_1n; a₂₁, a₂₂, a₂₃, ..., a_2n; ...; a_n1, a_n2, a_n3, ..., a_nn, viene espressa come prodotto tra la matrice L e la matrice U.

La matrice L ha elementi 1 sulla diagonale principale, elementi l₂₁, l₃₁, l₃₂, ..., l_n1, l_n2, l_n3, ..., e la matrice U ha elementi u₁₁, u₁₂, u₁₃, ..., u_1n; u₂₂, u₂₃, ..., u_2n; u₃₃, ..., u_3n; ...; u_nn.

Si può notare che la matrice L ha tutti uno sulla diagonale principale e che la prima riga della matrice U è uguale alla prima riga della matrice A.

Questa fattorizzazione è possibile se det(A_k) ≠ 0. Dove A_k è la sottomatrice di ordine k.

Gli elementi di L e U si ottengono dall’uguaglianza A = LU.

Come si procede con il calcolo?

Essendo che prima riga di U è sempre uguale alla prima riga di A. Per la seconda riga si avrà:

a₂₁ = l₂₁u₁₁ → l₂₁ = a₂₁ / u₁₁.

a₂₂ = l₂₁u₁₂ + l₂₂u₂₂ → u₂₂ = a₂₂ − l₂₁u₁₂.

a_2n = l₂₁u_1n + l₂₂u_2n → u_2n = a_2n − l₂₁u_1n.

Iterando si arriva alla conclusione che l’elemento a_mn è dato dalla somma dei prodotti degli elementi della riga m della matrice L per gli elementi della colonna n della matrice U, il normale prodotto riga per colonna.

In definitiva possiamo dire che:

u_ij = a_ij − Σ_k=1ⁿ⁻¹ l_iku_kj.

Quindi in generale si avrà che:

u_ij = a_ij − Σ_k=1ⁱ⁻¹ l_iku_kj, j = i, ..., n.

l_ij = (a_ij − Σ_k=1^j−1 l_iku_kj) / u_jj, j = 1, ..., i − 1.

4.3 Applicazioni della fattorizzazione

Soluzione di un sistema lineare.

Ax = b −→ LUx = b.

Ly = b, sistema triangolare inferiore.

Ux = y, sistema triangolare superiore.

Una volta fattorizzata A la soluzione del sistema si ottiene risolvendo i due sistemi triangolari.

Se si devono risolvere più sistemi lineari aventi la stessa matrice dei coefficienti, si fattorizza una sola volta A e per ogni vettore b si risolvono i due sistemi triangolari.

Calcolo del determinante di A.

det(A) = det(LU) = det(L) det(U), teorema di Binet.

Sapendo che det(L) = 1 si avrà che:

det(A) = det(LU) = det(U) = Π_k=1ⁿ u_kk = a₁₁⁽¹⁾ a₂₂⁽²⁾ ... a_nn⁽ⁿ⁾.

4.4 Applicazioni della fattorizzazione pt. 2

La matrice inversa di una matrice A non singolare è la matrice A⁻¹ tale che:

A A⁻¹ = I.

Si ricordi che I = E₁, E₂, ..., E_n, dove E_i sono vettori della base canonica.

Quindi:

A⁻¹ = X₁, X₂, ..., X_n.

A X_i = E_i.

Dove X_i sono le colonne di X che è la matrice inversa.

Conoscendo la fattorizzazione di A basta risolvere gli n sistemi lineari, cioè:

L Y_i = E_i, U X_i = Y_i, i = 1, 2, ..., n.

4.5 MEG e fattorizzazione LU

Il metodo di eliminazione gaussiana si basa sull’idea di ridurre il sistema Ax = b ad un sistema equivalente della forma Ux = b, dove U è triangolare superiore e b è un nuovo termine noto. Questo sistema può essere risolto con il metodo delle sostituzioni all’indietro.

Indichiamo il sistema originario come A⁽¹⁾x = b⁽¹⁾. Per ottenere il sistema equivalente utilizzeremo il metodo di Gauss.

Consideriamo una matrice A ∈ R, non singolare, e supponiamo che l’elemento diagonale a₁₁ sia non nullo.

Quindi:

m_i1 = a_i1⁽¹⁾ / a₁₁⁽¹⁾, i = 1, 2, ..., n.

Dove a_ij⁽¹⁾ sono gli elementi di A⁽¹⁾. È possibile eliminare l’incognita x₁ da tutte le righe successive alla prima sottraendo dalla riga i-esima con i = 2, 3, ..., n la prima riga moltiplicata per m_i1, eseguendo la stessa operazione per il termine noto. In generale si avrà che:

a_ij⁽²⁾ = a_ij⁽¹⁾ − m_i1a_1j⁽¹⁾, i = 2, 3, ..., n.

b_i⁽²⁾ = b_i⁽¹⁾ − m_i1b₁⁽¹⁾, i = 2, 3, ..., n.

Dove b_i⁽¹⁾ sono le componenti di b⁽¹⁾ si avrà un sistema:

a₁₁⁽¹⁾x₁ + a₁₂⁽¹⁾x₂ + ... + a_1n⁽¹⁾x_n = b₁⁽¹⁾.

0x₁ + a₂₂⁽²⁾x₂ + ... + a_2n⁽²⁾x_n = b₂⁽²⁾.

...

0x₁ + a_n2⁽²⁾x₂ + ... + a_nn⁽²⁾x_n = b_n⁽²⁾.

Che indicheremo con A⁽²⁾x = b⁽²⁾, equivalente a quello di partenza.

In modo analogo, possiamo trasformare il sistema in modo da eliminare l’incognita x_i dalla righe i + 1, ..., n.

In generale si otterrà:

A^(k)x = b^(k), 1 ≤ k ≤ n.

Assunto che a_ii⁽ⁱ⁾ ≠ 0 per i = 1, ..., k − 1. È evidente che per k = n si otterrà un sistema triangolare superiore.

A⁽ⁿ⁾x = b⁽ⁿ⁾.

a₁₁⁽¹⁾x₁ + a₁₂⁽¹⁾x₂ + ... + a_1n⁽¹⁾x_n = b₁⁽¹⁾.

0x₁ + a₂₂⁽²⁾x₂ + ... + a_2n⁽²⁾x_n = b₂⁽²⁾.

0x₁ + 0x₂ + ... + a_nn⁽ⁿ⁾x_n = b_n⁽ⁿ⁾.

Possiamo indicare con U la matrice triangolare superiore A⁽ⁿ⁾. Gli a_kk^(k) vengono detti elementi pivotiali e devono essere non nulli. Per evidenziare le formule che consentono di trasformare il sistema k-esimo in quello k+1-esimo, per k = 1, ..., n-1 si assume che a_kk^(k) ≠ 0 e definiamo il moltiplicatore:

m_ik = a_ik^(k) / a_kk^(k), i = k + 1, ..., n.

Quindi:

a_ij^(k+1) = a_ij^(k) − m_ika_kj^(k), i, j = k + 1, ..., n.

b_i^(k+1) = b_i^(k) − m_ikb_k^(k), i = k + 1, ..., n.

In generale:

M₁A⁽¹⁾ = A⁽²⁾.

M₂A⁽²⁾ = A⁽³⁾.

A⁽³⁾ = M₂A⁽²⁾ = M₂M₁A⁽¹⁾.

...

U = A⁽ⁿ⁾ = M_n−1M_n−2...M₁A⁽¹⁾.

Quindi si può scrivere:

U = A⁽ⁿ⁾ = (M_n−1M_n−2...M₁)A.

(M_n−1M_n−2...M₁)⁻¹ U = (M_n−1M_n−2...M₁)⁻¹ (M_n−1M_n−2...M₁)A⁽¹⁾.

(M_n−1M_n−2...M₁)⁻¹ U = A⁽¹⁾.

Ponendo (M_n−1M_n−2...M₁)⁻¹ = L si ottiene:

LU = A⁽¹⁾.

Esempio da pagina 17 Slides 4.

4.6 Matrice simmetrica definita positiva

Sia A ∈ M(n, R), una matrice di ordine n a coefficienti reali e sia x ∈ Rⁿ un generico vettore riga. La matrice è definita positiva se il prodotto x^TAx è un numero positivo per ogni x diverso dal vettore nullo.

A = A^T e x^TAx > 0, ∀x ∈ Rⁿ, x ≠ 0.

4.7 Criterio di Sylvester

Una matrice simmetrica è definita positiva se e solo se det(A_k) per k = 1,2,...,n e dove A_k sono le sottomatrici principali dell’ordine k.

4.8 Fattorizzazione di Cholesky

Data A = A^T definita positiva, si avrà che:

A = L L^T.

Per questo algoritmo basta fare il prodotto e uguagliare il risultato all’elemento a_ij.

Formule generali:

l_jj = √(a_jj − Σ_k=1^j−1 l_jk²).

l_ij = (a_ij − Σ_k=1^j−1 l_ikl_jk) / l_jj.

Esempio da pagina 22 Slides 4.

Codice pagina 25 Slides 4.

4.9 Matrici sparse

Le matrici sparse sono matrici con predominanza di elementi nulli. Non esiste una definizione rigorosa di sparsità, la precedente è solo qualitativa.

La sp