Calcolo Numerico

Appunti per il corso di calcolo numerico

Ingegneria informatica e delle telecomunicazioni

Anno accademico 2004-2005. Si precisa che i seguenti appunti sono estratti da quelli del corso di calcolo numerico della Professoressa R. Morandi. La presente versione ne costituisce infatti una versione rivista per i corsi di Ingegneria Informatica e delle Telecomunicazioni tenuti dalla Dott.ssa Costanza Conti. È importante che lo studente tenga presente che tali appunti costituiscono solo delle linee guida per la preparazione dell’esame. Tutti gli argomenti saranno discussi più in dettaglio a lezione (incluso gli esercizi relativi ad ogni argomento) e comunque sono rintracciabili in uno dei testi consigliati. Gli eventuali errori contenuti in questi appunti (di stampa e non) saranno senz’altro segnalati a lezione.

Testi consigliati

• D. Bini, M. Capovani, O. Menchi: Metodi Numerici per l’Algebra Lineare, Zanichelli, Bologna, 1988.
• D. Bini, R. Bevilacqua, M. Capovani, O. Menchi: Metodi Numerici, Zanichelli, Bologna, 1992.
• G. Monegato: Fondamenti di Calcolo Numerico, Levrotto e Bella, Torino, 1998.

Algoritmi

Un algoritmo è:

• Una successione finita di istruzioni assegnate in modo non ambiguo
• La cui esecuzione consenta di passare da una situazione iniziale (dati) ad una situazione finale (risultati)
• In un tempo finito

Requisiti fondamentali

• Generalità (adatto a risolvere una classe di problemi)
• Ottimalità (rispetto al tempo, al numero di operazioni, alla stabilità, ecc.)

Calcolare la somma di n numeri

1. Leggi: n, a₁, a₂, ..., a_n
2. Poni s = 0
3. Per i = 1, 2, ..., n, poni s = s + a_i
4. Scrivi s
5. Stop

Calcolare il prodotto di n numeri

1. Leggi: n, a₁, a₂, ..., a_n
2. Poni p = 1
3. Per i = 1, ..., n, poni p = p × a_i
4. Scrivi p
5. Stop

Determinare N, N > 0 mediante un procedimento iterativo

1. Individua m, n > 0 tali che 2^2m < N < n^m+n
2. Calcola c = 2
3. Se c = N, allora c = N → Stop
4. Se c < N, sostituire c a m e tornare a 2
5. Se c > N, sostituire c a n e tornare a 2

Nota Bene: il procedimento può non finire (es. N = 2), quindi non è un algoritmo

Criterio di arresto

Decidiamo di fermarci se 2|c - N| ≤ ε, ε > 0 è una tolleranza. C è un'approssimazione di N nei limiti del criterio di arresto scelto e della tolleranza scelta ε.

Algoritmo

1. Leggi: N, m, n, ε
2. Calcola c = (m + n)/2
3. Se |c - N| ≤ ε, vai a 5
4. Se c < N, poni m = c e vai a 2. Altrimenti poni n = c e vai a 2
5. Scrivi c
6. Stop

Nota Bene: rischio: [M, N] grande e/o ε piccolo potrebbero implicare troppo tempo di calcolo, quindi può essere una buona idea fissare un numero massimo di iterazioni.

Calcolo del valore di un polinomio in un punto α

P(x) = a₀ + a₁x + a₂x² + ... + a_nxⁿ

Algoritmo (1)

1. Leggi: n, a₀, a₁, ..., a_n; α
2. Poni p = 1
3. Poni s = a₀
4. Per i = 1, 2, ..., n, poni p = α · p, s = s + a_i · p
5. Scrivi s
6. Stop

Operazioni: 2n moltiplicazioni, n addizioni

Formulazione di Horner

P(x) = (...((a_nx + a_n-1)x + a_n-2)... + a₁)x + a₀

Algoritmo (2)

1. Leggi: n, a₀, a₁, ..., a_n; α
2. Poni s = a_n
3. Per i = n - 1, n - 2, ..., 0, poni s = s · α + a_i
4. Scrivi s
5. Stop

Operazioni: n moltiplicazioni, n addizioni

Rappresentazione dei numeri

Sistemi di numerazione

• Binario: cifre 0, 1
• Ottale: cifre 0, 1, 2, 3, 4, 5, 6, 7
• Decimale: cifre 0, 1, 2, 3, 4, 5, 6, 7, 8, 9

Notazione posizionale:

Es. (123)₁₀ = 1 × 10² + 2 × 10¹ + 3 × 10⁰

In un elaboratore i numeri sono rappresentati in una base diversa da β = 10, quindi in generale sono necessari algoritmi di conversione.

Rappresentazione dei numeri interi

(a_m, a_m-1, ..., a₁, a₀), con a_i, i = 1, ..., m cifre del sistema di numerazione 0 ≤ a_i < β.

Numeri reali

Un numero reale x si rappresenta nella forma di virgola mobile normalizzata (V.M.N.):

b^±0.a₁a₂...a_m β, a₁ ≠ 0

Esempi:

• 0.1471 → 0.1471 × 10²
• 337.1 → 0.371 × 10³
• 0.0371 → 0.371 × 10^-1

Nota: x = 0 ha per convenzione mantissa = 0, caratteristica = 0.

Rappresentazione in macchina dei numeri interi

Un elaboratore dispone di un certo numero di bit per rappresentare i numeri (interi e reali).

Numeri interi:

Locazioni di 16/32 bit (più comuni)

• 1 bit per il segno
• Avendo a disposizione t bit, solo t-1 bit sono per la rappresentazione vera e propria

Qual è il massimo numero N rappresentabile con t bit?

Es: β = 10; N = 947584

t = 7 → 0 9 4 7 5 8 4 | {z } t = 5 → overflow

N_max = 10^t-1 - 1

Rappresentazione in macchina dei numeri reali

Segno, caratteristica, mantissa. Numeri reali di macchina: numeri con mantisse e caratteristiche rappresentabili esattamente con i bit a disposizione.

Il limite sulla caratteristica delimita l’ordine di grandezza dei dati.

Il limite sulla mantissa caratterizza la precisione.

Precisione

Il numero di cifre che l’elaboratore riserva alla mantissa dei numeri reali (m) determina la precisione con la quale l’elaboratore lavora. Se x ≠ 0 è un numero reale tale che la sua rappresentazione in v.m.n. richiede più di m cifre di mantissa:

x = ±0.a₁a₂...a_ma_m+1...a_m+i β, a₁ ≠ 0

In questo caso x viene approssimato, cioè viene sostituito da un numero reale ad esso vicino che è chiamato fl(x) (floating di x).

La determinazione di fl(x) può avvenire mediante:

• Troncamento: fl(x) = ±0.a₁a₂...a_m β
• Arrotondamento: fl(x) = ±0.a₁a₂...(a_m + 1) β, se a_m+1 ≥ β/2

Esempio

β = 10, m = 5

x = 0.9345781 × 10³

Tronc. → fl(x) = 0.93457 × 10³

Arr. → fl(x) = 0.93458 × 10³

Errori nella rappresentazione

Si definisce:

• Errore assoluto: e_A = |x - fl(x)|
• Errore relativo: e_R = |x - fl(x)| / |x|

Esempio

β = 10, m = 4

Tronc/Arr

1. x = 0.5895 × 10⁵, fl(x) = 0.5895 × 10⁵
2. |x - fl(x)| = 0.00004 × 10⁵ = 0.4 × 10^-3
3. e_A = 0.4 × 10^-3
4. e_R = 0.4 × 10^-3 / 0.58954 × 10⁵

Il numero di cifre che l’elaboratore riserva alla mantissa dei numeri reali (m) determina la precisione con la quale l’elaboratore lavora. Se x ≠ 0 è un numero reale tale che la sua rappresentazione in v.m.n. richiede più di m cifre di mantissa:

x = ±0.a₁a₂...a_ma_m+1...a_m+i β, a₁ ≠ 0

In questo caso x viene approssimato, cioè viene sostituito da un numero reale ad esso vicino che è chiamato fl(x) (floating di x).

La determinazione di fl(x) può avvenire mediante:

• Troncamento: fl(x) = ±0.a₁a₂...a_m β
• Arrotondamento: fl(x) = ±0.a₁a₂...(a_m + 1) β, se a_m+1 ≥ β/2

Viene fatto lo stesso tipo di operazione su x e y con lo stesso risultato. L’errore relativo è lo stesso, ma l’errore assoluto non è lo stesso.

L’errore assoluto è influenzato dall’ordine di grandezza del dato. L’errore relativo non è influenzato dall’ordine di grandezza del dato. L’errore relativo dà indicazioni sull’approssimazione operata sulla mantissa del dato, ovvero sulla accuratezza (o precisione) con cui il dato è approssimato.

L’errore relativo non dipende dalla caratteristica, ma dipende dalla mantissa.

Precisione di macchina

Limitazione superiore degli errori relativi:

• Se si opera per troncamento: |x - fl(x)| / |x| ≤ β^1-m, ∀x ≠ 0
• Se si opera per arrotondamento: |x - fl(x)| / |x| ≤ β^1-m / 2, ∀x ≠ 0

La quantità β^1-m o β^1-m / 2 è detta precisione di macchina e viene di solito indicata con il simbolo ε_m. Più avanti vedremo un algoritmo per determinarla.

Aritmetica finita

Dati e risultati sono memorizzati con un numero finito di cifre (m) di mantissa e qualunque operazione viene effettuata con un numero finito di cifre (m) di mantissa.

Nota: risultati temporanei delle operazioni sono memorizzati in apposite locazioni con più cifre di mantissa, ma in ogni caso il numero è finito (si effettua un troncamento/arrotondamento).

Le 4 operazioni effettuate su numeri di macchina non producono necessariamente un numero di macchina, quindi il risultato deve essere trasformato nel suo floating.

Quindi:

• x + y si scrive x ⨁ y e fl(fl(x) + fl(y))
• x - y si scrive x ⨉ y e fl(fl(x) - fl(y))
• x × y si scrive x ⨀ y e fl(fl(x) × fl(y))
• x/y si scrive x ⟋ y e fl(fl(x) / fl(y))

Operazioni in macchina

Somma algebrica di due numeri reali

1. Si trasforma il numero con caratteristica minore in modo che i due numeri abbiano la stessa caratteristica (uno dei due perde la forma v.m.n.);
2. Si sommano le mantisse (lasciando invariate le caratteristiche);
3. Si ricava il floating del risultato (si rinormalizza il numero troncando o arrotondando se necessario)

Prodotto/Divisione di due numeri reali

1. Si producono/dividono le mantisse e si sommano/sottraggono le caratteristiche
2. Si fa il floating del risultato (si rinormalizza il numero troncando o arrotondando se necessario)

β = 10, m = 4, b = 2, tronc. x = 10^-1, y = -0.000540_^-3

fl(x) = 0.1000 × 100, fl(y) = -0.5400 × 100
x ⨁ y = 0.1000 × 100 ⨉ 0.0005 × 100 = 0.0995 × 100 = 0.9950 × 10-3
x ⨀ y = 0.1000 × 100 ⟋ (-0.0005) × 10-3 = -0.5 × 10-4

Alcune conseguenze della precisione finita

Per le operazioni di macchina valgono ancora le proprietà delle quattro operazioni aritmetiche? Spesso no.

• Non vale la proprietà associativa
• Non vale la proprietà distributiva
• Vale la proprietà commutativa

Formule o algoritmi matematicamente equivalenti (che portano allo stesso risultato se applicati in aritmetica esatta) possono produrre risultati diversi in aritmetica finita. Lo studio degli errori di arrotondamento e della loro propagazione attraverso gli algoritmi è di fondamentale importanza per poter interpretare e valutare i risultati di un qualunque algoritmo che operi con numeri reali.

Esempio

β = 10, m = 4, fl(x) per troncamento

ε = 10-3 ← ε = β-m ma = 2000, b = 2.5, c = 7.84
(a ⨁ b) ⨁ c = (0.2000 × 104 ⨁ 0.2500 × 104) ⨁ 0.7800 × 101
= (0.2000 × 104 ⨁ 0.00025) ⨁ 0.7800 × 101
= 0.2002 × 104 ⨁ 0.7800 × 101
= 0.2002 × 104 ⨁ 0.00078
a ⨁ (b ⨁ c) = 0.2000 × 104 ⨁ (0.2500 × 104 ⨁ 0.7800 × 101)
= 0.2000 × 104 ⨁ 0.1030 × 104
= 0.2000 × 104 ⨁ 0.00103 × 103
= 0.2010 × 104

(a + b) + c = a + (b + c) = 2010.3; (a ⨁ b) ⨁ c ≠ a ⨁ (b ⨁ c)

Note: errori assoluti 2010.3 - 2009 = 1.3; errori relativi 1.3/2010.3 ≈ 10^-3

Nota: risultati diversi ma accettabili nei limiti della precisione usata.

Importante

Il concetto di uguaglianza va modificato quando si lavora sui numeri reali in precisione finita. Risultati diversi possono essere considerati uguali nei limiti della precisione usata. Due numeri sono uguali quando hanno uguale caratteristica e uguale mantissa. Due numeri sono da considerarsi uguali quando la loro differenza è piccola rispetto alla precisione di macchina.

Ovvero due numeri a, b in precisione finita sono da considerarsi uguali quando |a - b| ≤ ε dove ε è un numero, detto tolleranza, il cui valore dipende dalla precisione di macchina.

Quindi

Nell’esempio del punto medio i due numeri a = 0.651 e b = 0.653 potevano essere considerati uguali. Infatti |b - a| = 2 × 10^-3 e ε = 10^-m.

Nota Bene: gli stessi numeri non potrebbero essere considerati uguali se fosse ε = 10^-8.

Se il risultato di un algoritmo in precisione infinita è zero e l’algoritmo realizzato in precisione finita dà un risultato diverso da zero ma piccolo, non possiamo concludere che l’elaboratore abbia sbagliato.

Se un algoritmo prevede un controllo sull’uguaglianza tra due reali dobbiamo fare molta attenzione nel realizzarlo in un programma. a = b → |a - b| ≤ ε

Il valore ε deve essere scelto tenendo conto della precisione di macchina e, in generale, da considerazioni indotte dal contesto.

Algoritmo per il calcolo di ε_m

La conoscenza di ε è fondamentale per valutare i risultati e scegliere le tolleranze. Serve un algoritmo per calcolarla. Idea su cui può essere basato l’algoritmo:

β = 10, m = 8, troncamento → ε = 10-7 = β-m
calcoliamo fl(1 + 10-6), fl(1 + 10-7), fl(1 + 10-8)
fl(1 + 10-6) = 0.10000000 × 101 ⨁ 0.10000000 × 101 = 0.10000010 × 101
fl(1 + 10-7) = 0.10000000 × 101 ⨁ 0.10000000 × 101 = 0.10000001 × 101
fl(1 + 10-8) = 0.10000000 × 101 ⨁ 0.10000000 × 101 = 0.100000001 × 101

ε è la più piccola potenza della base β = 10 che viene sentita in macchina se sommata a 1.

Algoritmo per determinare ε in base β

fl(1 + β^-p) ≠ fl(1) → ε = β^-p-1

1. Poni ε_m = 1
2. Poni ε_m = β × ε_m
3. Poni a = 1 + ε_m
4. Se a > 1 vai a 2. Altrimenti vai a 5.
5. Poni ε_m = β × ε_m
6. Stampa ε_m

Considerazioni sulle quattro operazioni aritmetiche

Supponiamo di effettuare operazioni esatte. Prodotto:

fl(x1) fl(x2) = x1x2(1 + ε1) × x2(1 + ε2)
= x1x2 - x1x2ε1 - x1x2ε2