NUMERI MACCHINA^{∗}

A. SOMMARIVA^{†}

Conoscenze richieste. Equazioni di secondo grado, successioni, formula di Taylor in una e due variabili, integrali di Riemann.

Conoscenze ottenute. In questa sezione analizziamo la rappresentazione floating point di un numero. Di seguito introduciamo il concetto di errore assoluto e relativo, condizionamento, stabilità e complessità computazionale con esempi in Matlab.

1. Numeri macchina. Fissato un numero naturale β > 1 è possibile vedere che ogni numero reale x ≠ 0 ha una unica rappresentazione del tipo

x = sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k}, 0 < d_{1}, 0 ≤ d_{k} ≤ β − 1, e ∈ ℤ

dove

sign(y) = \{ \begin{matrix}1, y > 0 \\ −1, y < 0 \\ 0, y = 0. \end{matrix}

è la funzione segno. La particolarità della sommatoria è ovviamente tale che x potrebbe essere non rappresentabile esattamente dal calcolatore se esistono infiniti d_{k} > 0, in quanto il calcolatore è capace di immagazzinare solo un numero finito di tali d_{k}. A tal proposito, definiamo il sottoinsieme F (β, t, L, U) di ℝ, costituito da quei numeri y che sono 0 oppure sono rappresentati con il sistema a virgola mobile normalizzata o floating point normalizzato da

y = sign(y) · (0.d_{1} … d_{t})β^{e} = sign(y) · β^{e} ∑_{k=1}^{t} d_{k}β^{−k}, d_{1} > 0, 0 ≤ d_{k} ≤ β − 1

con β > 1 numero naturale detto base, t numero prefissato di cifre di mantissa, e l’esponente intero e ∈ ℤ tale che L ≤ e ≤ U. Ogni elemento di F (β, t, L, U) è detto numero macchina. La parola normalizzato sottolinea che d_{1} > 0.

La base utilizzata dalla maggior parte dei computer attuali è β = 2, ma fino a non molti anni fa alcuni calcolatori usavano ancora la base β = 10 come l’HP-15C (≈ 1980) o β = 16 nel caso dell’IBM 3033 (≈ 1977) .

Di ogni numero macchina si immagazzinano una cifra per il segno, le cifre della mantissa e le cifre dell’esponente. Ognuna di queste cifre è detta bit. Ogni numero è registrato in vettore detto parola di un certo numero di M bits. Fissato M per la precisione singola, nel caso di precisione doppia questo è usualmente composto da 2 parole di M bits. Classici esempi sono per la precisione singola F (2, 24, −126, 128) in cui ogni numero macchina occupa 32 bit (o 4 byte essendo un byte pari a 8 bit) e per la doppia F (2, 53, −1022, 1024) in cui ogni numero macchina occupa 64 bit.

∗Ultima revisione: 25 settembre 2011 †Dipartimento di Matematica Pura ed Applicata, Università degli Studi di Padova, stanza 419, via Trieste 63,

35121 Padova, Italia (alvise@euler.math.unipd.it). Telefono: +39-049-8271350.

Osserviamo che per e ≥ t con t numero di cifre della mantissa, essendo β > 1 un numero naturale, allora

±β^{e} ∑_{k=1}^{t} d_{k}β^{−k} = ± ∑_{k=1}^{t} d_{k}β^{e−k}

è un numero intero, essendo e − k > 0. Di conseguenza i numeri macchina più grandi in valore assoluto sono interi. Si vede facilmente che tutti gli interi in valore assoluto inferiori o uguali a β^{U} · (1 − β^{−t}) sono in particolare nell’insieme F (β, t, L, U).

L’insieme dei numeri macchina F (β, t, L, U) ha un numero finito di elementi e si può dimostrare che la sua cardinalità è

2(β − 1)β^{t−1}(U − L + 1) + 1.

Se x è il numero reale

x = sign(x)(0.d_{1}, … d_{t}, …)β^{e} = sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k}, d_{1} > 0

allora con fl(x) denotiamo

fl(x) = sign(x)(0.d_{1}, …, d_{t−1}, \overline{d_{t}})β^{e} = sign(x)β^{e}(∑_{k=1}^{t−1} d_{k}β^{−k} + \overline{d_{t}}β^{−t})

in cui se si effettua il troncamento

\overline{d_{t}} = d_{t}, d_{1} > 0

mentre se viene effettuato l’arrotondamento.

\overline{d_{t}} = \{ \begin{matrix}d_{t}, d_{t+1} < β/2 \\ d_{t} + 1, d_{t+1} ≥ β/2 \end{matrix}

La scelta tra troncamento o arrotondamento è eseguita a priori dal sistema numerico utilizzato dal calcolatore. Attualmente la scelta di arrotondare è la più diffusa.

Se l’esponente e del numero x è minore di L si commette un’errore di underflow, mentre se è maggiore di U si commette un’errore di overflow. Se invece e ∈ [L, U], x = ±β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k} e per qualche t^{∗} > t si ha d_{t^{∗}} > 0 allora si commette un errore di troncamento o arrotondamento a seconda si effettui l’uno o l’altro.

Essendo un numero macchina del tipo

x = ±β^{e} ∑_{k=1}^{t} d_{k}β^{−k}

con d_{1} > 0 e 0 ≤ d_{k} ≤ β − 1 si vede subito che

|x| = β^{e} ∑_{k=1}^{t} d_{k}β^{−k} ≤ β^{e} ∑_{k=1}^{t}(β − 1)β^{−k}

= β^{e}(β − 1) ∑_{k=1}^{t} β^{−k} = β^{e}(β − 1) (\frac{1 − β^{−(t+1)}}{1 − β^{−1}} − 1)

= β^{e}(β − 1) (\frac{β(1 − β^{−(t+1)})}{β − 1} − 1)
= β^{e}(β − 1) (\frac{β(1 − β^{−(t+1)}) − (β − 1)}{β − 1})
= β^{e}(1 − β^{−t}) < β^{e}  (1.1)

In altri termini è minore di 1 la somma relativa alla mantissa, cioè ∑_{k=1}^{t} d_{k}β^{−k} < 1. Poichè , dai calcoli eseguiti si vede che β^{e}(1 − β^{−t}) è un numero macchina, si deduce

che il massimo numero macchina rappresentabile è quindi

β^{U}(1 − β^{−t}).

Si vede facilmente che il più piccolo numero macchina (e quindi normalizzato!) positivo rappresentabile è β^{L−1}.

Si dice precisione di macchina, quel numero ξ_{m} che rappresenta la distanza tra 1 ed il successivo numero in virgola mobile. Calcoliamo tale valore. Necessariamente il numero successivo a 1 ha per mantissa d_{1} = 0, …, d_{t−1} = 0, d_{t} = 1 e esponente e = 1 (perchè deriva da (1.1)?) e quindi è x = β^{1} · β^{−t} = β^{1−t}. Si sottolinea che la precisione di macchina non coincide in generale con il più piccolo numero rappresentabile dal calcolatore. Più in generale la distanza tra x = sign(x)β^{e} ∑_{k=1}^{t} d_{k}β^{−k} e il floating point successivo è uguale a sign(x)β^{e−t} e quindi i numeri floating point non sono equispaziati.

Fissato un vettore da approssimare

x^{∗} = (x_{1}^{∗}, …, x_{n}^{∗}) ∈ ℝ^{n},

si definisce

1. errore assoluto tra x e x^{∗} la quantità

∥x − x^{∗}∥

dove

∥y∥ = \sqrt{∑_{i=1}^{n} y_{i}^{2}}, y = (y_{1}, …, y_{n});

2. errore relativo tra x e x^{∗} ≠ 0 la quantità

\frac{∥x − x^{∗}∥}{∥x^{∗}∥}.

Si noti che se ∥x^{∗}∥ = 0, cioè x^{∗} = 0 per la proprietà delle norme [1, p.481], allora non ha senso la scrittura propria dell’errore relativo. Inoltre se il vettore ha un solo componente, allora la norma ∥ · ∥ coincide con l’usuale valore assoluto | · |.

Sia

x = sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k} ≠ 0, d_{1} > 0

un numero reale e

fl(x) = sign(x)(0.d_{1}, …, d_{t−1}, \overline{d_{t}})β^{e} = sign(x)β^{e}(∑_{k=1}^{t−1} d_{k}β^{−k} + \overline{d_{t}}β^{−t}), d_{1} > 0

la sua approssimazione in virgola mobile.

L’errore relativo compiuto nell’approssimare x con fl(x), visto che |x| ≥ |d_{1}β^{e−1}| ≥ |β^{e−1}| poichè d_{1} > 0, verifica

\frac{|x − fl(x)|}{|x|} ≤ \frac{|sign(x)β^{e} ∑_{k=1}^{+∞} d_{k}β^{−k} − (sign(x)β^{e}(∑_{k=1}^{t−1} d_{k}β^{−k} + \overline{d_{t}}β^{−t}))|}{|β^{e−1}|}

= β · |∑_{k=t}^{+∞} d_{k}β^{−k} − \overline{d_{t}}β^{−t}|  (1.2)

Nel caso si effettui il troncamento abbiamo d_{t} = \overline{d_{t}} e quindi, essendo ξ_{M} = β^{1−t} la precisione di macchina

β · |∑_{k=t}^{+∞} d_{k}β^{−k} − \overline{d_{t}}β^{−t}| = β · | ∑_{k=t+1}^{+∞} d_{k}β^{−k}| ≤ β · (β · β^{−t−1}) = β^{1−t} = ξ_{M}

mentre se si effettua l’arrotondamento abbiamo dopo alcuni calcoli che

\frac{|x − fl(x)|}{|x|} < \frac{β^{1−t}}{2} = \frac{ξ_{M}}{2} .

Il minor errore relativo compiuto nell’approssimare x con fl(x) via arrotondamento, suggerisce perchè quest’ultimo sistema sia più utilizzato. Da questo si deduce che i numeri in virgola mobile non sono dunque un insieme denso, come sono invece i numeri reali, bensì un insieme discreto di valori sull’asse reale, in posizioni non equispaziate.

2. Facoltativo: Rappresentazione IEEE 754 dei numeri al calcolatore. La rappresentazione numerica nei moderni calcolatori rispetta in genere lo standard IEEE 754 che andiamo brevemente a descrivere.

Si intende per bit una cifra nel sistema (usualmente binario) usato [2, p. 48], [9, p. 47]. Nello standard IEEE 754, i numeri in virgola mobile a precisione singola sono memo-

rizzati in un vettore detto parola di 32 bits, mentre quelli a precisione doppia occupano due parole consecutive da 32 bits. Un numero non nullo-normalizzato si scrive come

x = (−1)^{S} · 2^{E−Bias} · (1 + F)  (2.1)

dove S = 0 oppure S = 1 determina il segno, E − Bias è l’esponente cui va sottratto il valore del BIAS ed F la mantissa con F = ∑_{k=1}^{+∞} a_{k}2^{−k}, con a_{k} = 0 oppure a_{k} = 1.

Evidentemente, in un calcolatore non si possono immagazzinare un numero infinito di a_{k}, e ci si limita ai cosidetti numeri macchina

x = (−1)^{S} · 2^{E−Bias} · (1 + F)  (2.2)

con E intero tale che |E| < 2^{N_{E}} per un certo prefissato N_{E} e F = ∑_{k=1}^{N_{F}} a_{k}2^{−k}, Alcuni esempi

Precisione N_{E} N_{F} Bias
Singola  8 23 127
Doppia  11 52 1023

Partiamo con un esempio semplice. Il numero 12 corrisponde in singola precisione a

[0 | 10000010 | 10000000000000000000000]

Cerchiamo di comprenderne il significato:

1. 0 determina il segno. Poichè (−1)^{0} = 1 il segno è +.

2. 10000010 determina il valore della potenza di 2. Essendo il BIAS uguale a 127 e (si legga da destra a sinistra!) E = 10000010 = 0 · 2^{0} + 1 · 2^{1} + 0 · 2^{2} + 0 · 2^{3} + 0 · 2^{4} + 0 · 2^{5} + 0 · 2^{6} + 1 · 2^{7} = 2 + 128 = 130, l’esponente vale 130 − 127 = 3. Quindi 2^{E−Bias} = 2^{3} = 8.

3. F = 10000000000000000000000 è la componente relativa alla mantissa. La quantità , letta da sinistra a destra, riguarda la parte decimale F in (2.2). Quindi F = 1 · 2^{−1} + 0 · 2^{−2} + … = 0.5. Ora (1 + F) si legge come 1 + F = 1.5.

In altre parole la rappresentazione di 12 è esatta ed il numero è rappresentato come 12 = 8 · 1.5.

Non convinti, proviamo un altro esempio, ad esempio 126 che corrisponde in singola precisione a

[0 | 10000101 | 11111000000000000000000]

1. 0 determina il segno. Poichè (−1)^{0} = 1 il segno è +.

2. 10000101 determina il valore. Essendo il BIAS uguale a 127 e (si legga da destra a sinistra!) E = 10000101 = 1 · 2^{0} + 0 · 2^{1} + 1 · 2^{2} + 0 · 2^{3} + 0 · 2^{4} + 0 · 2^{5} + 0 · 2^{6} + 1 · 2^{7} = 1 + 4 + 128 = 133, l’esponente vale 133 − 127 = 6. Quindi 2^{E−Bias} = 2^{6} = 64. Viene da chiedersi perchè introdurre il BIAS. La risposta è che rappresenta un alternativa all’introduzione di un bit di segno per l’esponente.

3. F = 11111000000000000000000 è la componente relativa alla mantissa. La quantità , letta da sinistra a destra, riguarda la parte decimale F in (2.2). Quindi F = 1 · 2^{−1} + 1 · 2^{−2} + 1 · 2^{−3} + 1 · 2^{−4} + 1 · 2^{−1} + … = 0.5 + 0.25 + 0.125 + 0.0625 + 0.03125 = 0.96875. Ora (1 + F) si legge come 1 + F = 1.96875.

In altre parole la rappresentazione di 126 è esatta ed il numero è rappresentato come 126 = 64 · 1.96875.

La distanza tra un numero in virgola mobile x^{−} ed il successivo x^{+}, è

2^{−N_{F}} · 2^{E−bias}.

Per convincersene, supponiamo che siano

x^{−} = (−1)^{S} · 2^{E−Bias} · (1 + F^{−}), x^{+} = (−1)^{S} · 2^{E−Bias} · (1 + F^{+}).

Allora, essendo F^{+} − F − = 2^{−N_{F}} (pensarci su!) facilmente

|x^{+} − x^{−}| = 2^{E−Bias} · (F^{+} − F −) = 2^{−N_{F}} · 2^{E−bias}.

La cardinalità dei numeri macchina in IEEE 754 è piuttosto elevata. Ad esempio, in doppia precisione sono

2^{64} ≈ 1.844674407370955e + 019.

La precisione di macchina, cioè quel numero che rappresenta la distanza tra 1 ed il successivo numero in virgola mobile. Per avere un’idea del suo ordine di grandezza, tale valore eps in Matlab 7.6.0 vale

>> format long
>> eps
ans =
      2.220446049250313e-16
>>

3. Facoltativo: Su eps, realmax, realmin, overflow e underflow. I numeri nella rappresentazione floating-point stanno in un intervallo limitato, quindi l’infinitamente piccolo e l’infinitamente grande non sono rappresentabili come numeri in virgola mobile, nel calcolatore. Quando si raggiunge un valore talmente piccolo da non essere pi distinto dallo zero, si parla di underflow, mentre quando si eccede il massimo numero rappresentabile, si parla di overflow. Il limite inferiore (che non sia denormalizzato!) è pari a 2^{1−bias} ed in Matlab è espresso dalla variabile realmin. Il limite superiore è pari a 2^{2_{E}^{N} −2−bias} ·(1+(1−2-Fbits)) ed in Matlab è espresso dalla variabile realmax. Verifichiamolo.

1. Dalla tabella, in precisione doppia, si ha che bias=1023 e quindi

2^{1−1023} = 2^{−1022} = 2.225073858507201e − 308

Se digitiamo nella shell di Matlab il comando realmin otteniamo proprio questo valore.

2. Dalla tabella, in precisione doppia, si ha che bias= 1023, N_{E} = 11, ed N_{F} = 52. Conseguentemente

realmax = 2^{2^{N_{E}} −2−bias} · (1 + (1 − 2^{N_{F}}))
= 2^{2^{11}−2−1023} · (1 + (1 − 2^{−52}))
= 2^{2^{11}−2−1023} · (2 − 2^{−52}) =
= 1.797693134862316e + 308

Per convincercene, digitiamo nella shell di Matlab/Octave

>> s=2^11-1025;
>> m=2-2^(-52)
m =
     2.0000
>> m=2-2^(-52);
>> t=(2^s)*m
t =
  1.7977e+308
>> realmax
ans =
  1.7977e+308
>> t-realmax
ans =
      0
>>

Quindi digitando realmax nella shell di Matlab otteniamo proprio questo valore. Osserviamo che in realtà realmin non sembra essere il numero più piccolo in valore

assoluto. Testiamolo sulla shell di Matlab.

>> help realmax
REALMAX Largest positive floating point number.
x = realmax is the largest floating point number representable
on this computer.  Anything larger overflows.
See also EPS, REALMIN.
Overloaded methods
help quantizer/realmax.m
>> a=realmax \cdot 2
a =
Inf
>> help realmin
REALMIN Smallest positive floating point number.
x = realmin is the smallest positive normalized floating point
number on this computer.  Anything smaller underflows or is
an IEEE "denormal".
See also EPS, REALMAX.
Overloaded methods
help quantizer/realmin.m
>> realmin/(2^51) == 0
ans =
0
>> realmin/(2^52) == 0
ans =
0
>> realmin/(2^53) == 0

ans =
      1
>>

Siccome in un costrutto logico il valore 1 corrisponde a true, si ricava che

realmin/(2^{52}) ≈ 4.940656458412465e − 324

è ancora un numero macchina, ma che fa parte dei cosidetti numeri denormalizzati (cf. [9], p.49) e cioè del tipo

x = (−1)^{S} · 2^{E−Bias} · (0.F).  (3.1)

Per curiosità vediamo il codice Matlab che implementa realmin:

function xmin = realmin
%REALMIN Smallest positive floating point number.
%   x = realmin is the smallest positive normalized floating point
%   number on this computer.  Anything smaller underflows or is
%   an IEEE "denormal".
%
%   See also EPS, REALMAX.
%   C. Moler, 7-26-91, 6-10-92.
%   Copyright 1984-2001 The MathWorks, Inc.
%   $Revision: 5.8 $  $Date: 2001/04/15 12:02:27 $
minexp = -1022;
% pow2(f,e) is f \cdot 2^e, computed by adding e to the exponent of f.
xmin = pow2(1,minexp);

Per quanto concerne pow2:

>> help pow2
 POW2   Base 2 power and scale floating point number.
     X = POW2(Y) for each element of Y is 2 raised to the power Y.
     X = POW2(F,E) for each element of the real array F and a integer
     array E computes X = F . \cdot  (2 .^E).  The result is computed
     quickly by simply adding E to the floating point exponent of E.
     This corresponds to the ANSI C function ldexp() and the IEEE
     floating point standard function scalbn().
     See also LOG2, REALMAX, REALMIN.
>>

4. Operazioni con i numeri macchina. Per quanto riguarda le operazioni aritmetiche fondamentali si può dimostrare che la somma, la divisione e la moltiplicazione producono un errore relativo piccolo, mentre la sottrazione può anche produrre un errore relativo grande rispetto al risultato (ciò avviene quando i due operandi sono molto vicini tra di loro e si ha dunque una perdita di cifre significative nel risultato). Più precisamente se fl(x) è il numero macchina che corrisponde a x, denotiamo con

x ⊕ y = fl(fl(x) + fl(y))  (4.1)
x ⊖ y = fl(fl(x) − fl(y))  (4.2)
x ⊗ y = fl(fl(x) · fl(y))  (4.3)
x ⊘ y = fl(fl(x) : fl(y))  (4.4)

e per ⊙ una delle operazioni precedentemente introdotte, cioè una tra ⊕, ⊖, ⊗, ⊘ corrispondenti a op cioè una tra +, −, ·, :, posto

ϵ_{x} = \frac{|x − fl(x)|}{|x|}  (4.5)
ϵ_{x,y}^{⊙} = \frac{|(x opy) − (x ⊙ y)|}{|x opy|}  (4.6)

Per prima cosa osserviamo che alcune proprietà caratteristiche di +, −, · e : non valgono per i corrispettivi ⊕, ⊖, ⊗, ⊘. Infatti per la somma ⊕ e il prodotto ⊗ vale la proprietà commutativa ma non vale quella associativa, e nemmeno quella distributiva.

Vediamo alcuni esempi. Mostriamo come non valga la proprietà associativa della somma. Consideriamo F (10, 4, −10, +10) con fl basato sul troncamento. Sia a = 2000, b = 2.5, c = 7.8. Allora:

(a ⊕ b) ⊕ c = 0.2002 · 10^{4} ⊕ 0.7800 · 10^{1} = 0.2009 · 10^{4}
a ⊕ (b ⊕ c) = 0.2000 · 10^{4} ⊕ 0.1030 · 10^{2} = 0.2010 · 10^{4}

e quindi (a ⊕ b) ⊕ c ≠ a ⊕ (b ⊕ c). Per quanto riguarda la proprietà distributiva, in F (10, 4, −10, +10), posti a = 0.800000 · 10^{9}, b = 0.500009 · 10^{4}, c = 0.500008 · 10^{4}

a ⊗ (b ⊖ c) = 0.800000 · 10^{9} ⊗ 0.1 · 10^{−2} = 0.800000 · 10^{9}
(a ⊗ b) ⊖ (a ⊗ c) = (0.800000 · 10^{9} ⊗ b = 0.500009 · 10^{4}) ⊖
⊖(0.800000 · 10^{9} ⊗ b = 0.500008 · 10^{4}) → overflow (4.7)

(l’overflow è dovuto alla sequenza di operazioni da compiere).

4.0.1. Errori nelle operazioni e loro propagazione. Consideriamo ora come si propagano gli errori in ⊕, ⊖, ⊗, ⊘. Dopo qualche tedioso conto (cf. [2, p.78])

ϵ_{x,y}^{⊕} ≤ |\frac{x}{x + y} |ϵ_{x} + |\frac{y}{x + y} |ϵ_{y}  (4.8)
ϵ_{x,y}^{⊖} ≤ |\frac{x}{x − y} |ϵ_{x} + |\frac{y}{x − y} |ϵ_{y}  (4.9)
ϵ_{x,y}^{⊗} ⪅ ϵ_{x} + ϵ_{y}  (4.10)
ϵ_{x,y}^{⊘} ≤ |ϵ_{x} − ϵ_{y}|  (4.11)

il che mostra il pericolo di perdita di accuratezza nella somma e nella sottrazione qualora rispettivamente x + y ≈ 0 e x − y ≈ 0 (fenomeno di cancellazione). Proviamo di seguito il caso della somma e della moltiplicazione, lasciando al lettore quelli di sottrazione e divisione. Con a ⪅ b si intende a inferiore o uguale circa a b.

Per quanto riguardo la somma, supposti x + y ≠ 0, x ≠ 0, y ≠ 0, e fl(fl(x) + fl(y)) = fl(x) + fl(y) (cioè fl(x) + fl(y) numero macchina, richiesta non eccessiva), abbiamo

ϵ_{x,y}^{⊕} = \frac{|(x + y) − (x ⊕ y)|}{|x + y|} = \frac{|(x + y) − fl(fl(x) + fl(y))|}{|x + y|} = \frac{|(x − fl(x)) + (y − fl(y))|}{|x + y|}
≤ \frac{|x|}{|x + y|} · \frac{|x − fl(x)|}{|x|} + \frac{|y|}{|x + y|} · \frac{|y − fl(y)|}{|y|} = |\frac{x}{x + y} |ϵ_{x} + |\frac{y}{x + y} |ϵ_{y}  (4.12)

Nel caso del prodotto, supposti x + y ≠ 0, x ≠ 0, y ≠ 0, e fl(fl(x) · fl(y)) = fl(x) · fl(y) (cioè fl(x) · fl(y) numero macchina), abbiamo

ϵ_{x,y}^{⊗} = \frac{|x · y − x ⊗ y|}{|x · y|} = \frac{|x · y − fl(fl(x) · fl(y))|}{|x · y|} = \frac{|x · y − fl(x) · fl(y)|}{|x · y|}
= \frac{|x · y − x · fl(y) + x · fl(y) − fl(x) · fl(y)|}{|x · y|}
≤ \frac{|x · (y − fl(y))| + |fl(y) · (x − fl(x))|}{|x · y|}
= \frac{|y − fl(y)|}{|y|} + \frac{|fl(y) · (x − fl(x))|}{|x · y|}
≈ ϵ_{x} + ϵ_{y}  (4.13)

dove si ricorda che

\frac{|fl(y)|}{|x · y|} ≈ \frac{1}{|x|}.

Nota 4.1. Osserviamo che dalla formula di Taylor in due variabili, posto per brevità di notazione f′_{x} = ∂f/∂x, f′_{y} = ∂f/∂y,

f(x, y) ≈ f(\overline{x}, \overline{y}) + f′_{x}(\overline{x}, \overline{y})(x − \overline{x}) + f′_{y}(\overline{x}, \overline{y})(y − \overline{y})  (4.14)

necessariamente, se x, y, f(x, y) ≠ 0 allora posti ϵ_{x} = |x − \overline{x}|/|x|, ϵ_{y} = |y − \overline{y}|/|y| gli errori relativi sui dati

\frac{|f(x, y) − f(\overline{x}, \overline{y})|}{|f(x, y)|} ⪅ \frac{|f′_{x}(\overline{x}, \overline{y})||x − \overline{x}||x|}{|f(x, y)||x|} + \frac{|f′_{y}(\overline{x}, \overline{y})||y − \overline{y}||y|}{|f(x, y)||y|}
= \frac{|f′_{x}(\overline{x}, \overline{y})||x|}{|f(x, y)|} · ϵ_{x} + \frac{|f′_{y}(\overline{x}, \overline{y})||y|}{|f(x, y)|} · ϵ_{y}  (4.15)

Da (4.15), si possono ottenere risultati simili ai precedenti per ⊕, ⊗.

5. Valutazione di una funzione. Nel valutare una funzione continua f : Ω ⊆ ℝ → ℝ in un punto x ∈ Ω bisogna tener conto di vari dettagli. Per prima cosa il dato x può essere per varie ragioni affetto da errori, tipicamente di misura o di rappresentazione al calcolatore e conseguentemente approssimato con un certo \overline{x}. Di conseguenza, spesso invece di f(x) si valuta f(\overline{x}). Ci si domanda se ci siano considerevoli differenze nel valutare la funzione in \overline{x} invece che in x.

Una funzione f risulta difficile da valutare al calcolatore nel punto x ≠ 0 in cui f(x) ≠ 0 qualora a piccoli errori (relativi) sui dati

|x − x_{c}|/|x|

corrispondano grandi errori (relativi) sui risultati

|f(x) − f(x_{c})|/|f(x)|.

Di conseguenza è importante discutere la quantità

K(f, x, x_{c}) = \frac{|f(x) − f(x_{c})|/|f(x)|}{|x − x_{c}|/|x|}.

La sua valutazione in un punto x si dice bencondizionata se a piccole variazioni relative sui dati corrispondono piccole variazioni sui risultati, malcondizionata se ciò non avviene. Se f è derivabile con continuità nel più piccolo intervallo I aperto contenente x ed x_{c}, per il teorema della media, essendo

f(x) − f(x_{c}) = f′(ξ) · (x − x_{c}), ξ ∈ I,

e quindi da |f(x) − f(x_{c})|/|x − x_{c}| = f′(ξ) ≈ f′(x) ricaviamo facilmente che

K(f, x, x_{c}) ≈ K(f, x) := \frac{|x · f′(x)|}{|f(x)|}.

La quantità K(f, x) si chiama condizionamento di f nel punto x. Più piccola risulta K(f, x) e meno la funzione amplifica un piccolo errore sul dato x. Naturalmente la funzione può essere bencondizionata su un certo dato x_{1} ∈ Ω e malcondizionata su un certo altro x_{2} ∈ Ω. Osserviamo che il condizionamento di una funzione non dipende dall’algoritmo con cui viene valutata, ma è inerente alla funzione stessa f e al dato x.

5.1. Esercizio sul condizionamento. Date le funzioni

f_{1}(x) = 1 − \sqrt{1 − x^{2}}
f_{2}(x) = 1 − x

si calcoli analiticamente il condizionamento

K(f_{1}, x) = \frac{|x · f′_{1}(x)|}{|f_{1}(x)|}

K(f_{2}, x) = \frac{|x · f′_{2}(x)|}{|f_{2}(x)|}

Utilizzando Matlab/Octave, si plottino in scala semilogaritmica i grafici di K(f_{1}, ·) e K(f_{2}, ·) nel set di punti x = −1 + 10^{(−3)}, −1 + 2 · 10^{(−3)}, …, 1 − 2 · 10^{(−3)}, 1 − 10^{(−3)},. Si può dire da questi dove la funzione è malcondizionata (cioè assume valori di condizionamento alti)?

6. Stabilità . Nella sezione precedente abbiamo analizzato il condizionamento di una funzione continua f in un punto x del suo dominio. In particolare, abbiamo supposto che f sia valutabile esattamente e quindi il bencondizionamento dipende esclusivamente da f e dal dato x e non dall’algoritmo utilizzato. Un algoritmo per la risoluzione di un certo problema si dice stabile se amplifica poco gli errori di arrotondamento introdotti dalle singole operazioni. Il problema non è necessariamente la valutazione di una funzione continua in un punto, ma può essere la risoluzione di una equazione di secondo grado, il calcolo di π greco, la valutazione di una successione, etc. Osserviamo inoltre che nel caso della valutazione in un punto x di una funzione f bencondizionata, si possono ottenere risultati non sufficientemente corretti se il problema è affrontato con un algoritmo instabile.

6.1. Consistenza, stabilità e convergenza. Vediamo i dettagli del concetto di stabilità e li legheremo ad altre due proprietà , quelle di consistenza e convergenza (cf.[9, p.39-41]). Supponiamo di voler risolvere il problema

F (x, d) = 0,  (6.1)

e che questo sia ben posto, cioè a piccole variazioni di d corrispondano piccole variazioni di x. Nella precedente scrittura, la variabile d è fissata e si cerca x che risolva il problema (6.1).

Un metodo numerico per la soluzione approssimata di F (x, d) costruisce una successione di problemi

F_{n}(x_{n}, d_{n}) = 0, n ≥ 1

dipendenti da n, col desiderio che sia

x_{n} → x, n → ∞.

Se il dato d è ammissibile per F_{n}, cioè ha senso risolvere F_{n}(x, d) = 0, si dice che il metodo è consistente se per x soluzione del problema F (x, d) = 0 si ha

F_{n}(x, d) = F_{n}(x, d) − F (x, d) → 0, n → ∞.

Il metodo si dice fortemente consistente se in particolare

F_{n}(x, d) = F_{n}(x, d) − F (x, d) = 0, ∀n.

Come detto in precedenza, un metodo è stabile (talvolta si dice ben posto) se per ogni n fissato esiste e sia unica la soluzione x_{n} del problema F (x_{n}, d_{n}) = 0 e che x_{n} dipenda con continuità dai dati cioè

∀η > 0, ∃K_{n}(η, d_{n}) : ∥δd_{n}∥ < η ⇒ ∥δx_{n}∥ ≤ K_{n}(η, d_{n})∥δd_{n}∥.

Il proposito di un metodo è ovviamente di costruire, attraverso problemi numerici

F_{n}(x_{n}, d_{n}) = 0,

una successione di x_{n} che si avvicinino sempre più a x al crescere di n, cioè il metodo sia convergente. In termini più rigorosi, se F (x(d), d) = 0 e F_{n}(x_{n}(d + δd_{n}), d + δd_{n}) = 0

∀ϵ > 0 ∃n_{0}(ϵ), ∃δ(n_{0}, ϵ) > 0 : ∀n > n_{0}(ϵ), ∀∥δd_{n}∥ < δ(n_{0}, ϵ)
⇒ ∥x(d) − x_{n}(d + δd_{n})∥ ≤ ϵ.  (6.2)

Il teorema di Lax-Richtmyer mostra come queste tre proprietà siano interconnesse, cioè

Teorema 6.1. Per un metodo numerico consistente, la stabilità è equivalente alla convergenza.

6.2. Esempio 1: Calcolo di una radice in una equazione di secondo grado. Vediamo un esempio concreto in cui l’introdurre una sottrazione potenzialmente pericolosa conduce effettivamente a problemi di instabilità della soluzione e come rimediare. Dato il polinomio di secondo grado x^{2} + 2 px − q, con \sqrt{p^{2} + q} ≥ 0 calcoliamo la radice

y = −p + \sqrt{p^{2} + q}.  (6.3)

Osserviamo che essendo \sqrt{p^{2} + q} ≥ 0 le radici

y = −p ± \sqrt{p^{2} + q}.  (6.4)

dell’equazione sono reali. La soluzione descritta in (6.3) è la maggiore delle 2, ed è non negativa se e solo se q ≥ 0. Consideriamo la funzione che fissato q ≥ 0, mappa un numero reale p in −p ± \sqrt{p^{2} + q}. Si osserva subito che (6.3) è potenzialmente instabile per p ≫ q a causa della sottrazione tra p e \sqrt{p^{2} + q}. A tal proposito, dopo averla implementata in Matlab, verificheremo numericamente la perdita di accuratezza per opportune scelte dei coefficienti p e q. Ripetiamo poi lo stesso tipo di indagine con una formula alternativa (e stabile) che si ottiene razionalizzando la formula (6.3). In altri termini

y = −p + \sqrt{p^{2} + q} = \frac{(−p + \sqrt{p^{2} + q})(p + \sqrt{p^{2} + q})}{(p + \sqrt{p^{2} + q})} = \frac{q}{(p + \sqrt{p^{2} + q})}  (6.5)

Ricordiamo ora che un problema si dice bencondizionato (o malcondizionato) a seconda che nel particolare contesto le perturbazioni sui dati non influenzino (o influenzino) eccessivamente i risultati. Nel caso di un algoritmo, per indicare un simile comportamento rispetto alla propagazione degli errori dovute alle perturbazioni sui dati, si parla di algoritmo bencondizionato (o algoritmo malcondizionato) anche se è più usuale il termine di stabilità [2, p. 66].

Seguendo [2, p. 10], [2, p. 78], il problema (e non l’algoritmo!) è bencondizionato per q > 0 e malcondizionato per q ≈ −p^{2}. Usando dei classici ragionamenti dell’analisi numerica si mostra che (cf. [10], p. 21, [3], p. 11)

1. il primo algoritmo (6.4) non è numericamente stabile qualora p ≫ q > 0;

2. il secondo algoritmo (6.5) è numericamente stabile qualora p ≫ q > 0. In [10, p.22], si suggerisce un test interessante per

p = 1000, q = 0.018000000081

la cui soluzione esatta è 0.9 · 10^{−5}. Secondo [3], p. 11 è notevole l’esperimento in cui

p = 4.999999999995 10^{+4}, q = 10^{−2}

avente soluzione esatta 10^{−7}. Si osservi che in entrambi i casi effettivemente p ≫ q. Passiamo ora all’implementazione e verifica numerica di questi due tests, osservando che i problemi relativi sono comunque ben condizionati.

Scriviamo un programma radicesecgrado.m in Matlab che illustri i due algoritmi.

% p=4.999999999995 \cdot 10^(+4); q=10^(-2); sol=10^(-7);
p=1000; q=0.018000000081; sol=0.9*10^(-5);
% ALGORITMO 1
s=p^2;
t=s+q;
if t >=0
   u=sqrt(t);
else
   fprintf('\n \t [RADICI COMPLESSE]');

end
s1=-p+u;
% ALGORITMO 2
s=p^2;
t=s+q;
if t >=0
   u=sqrt(t);
else
   fprintf('\n \t [RADICI COMPLESSE]');
end
v=p+u;
t1=q/v;
fprintf('\n \t [ALG.1] [1]: %10.19f',s1);
fprintf('\n \t [ALG.2] [1]: %10.19f',t1);
if length(sol) > 0 & (sol ~= 0)
     relerrs1=abs(s1-sol)/abs(sol);
     relerrt1=abs(t1-sol)/abs(sol);
     fprintf('\n \t [REL.ERR.][ALG.1]: %2.2e',relerrs1);
     fprintf('\n \t [REL.ERR.][ALG.2]: %2.2e',relerrt1);
end

Digitiamo quindi da shell Matlab/Octave il comando radicesecgrado e otteniamo

>> radicesecgrado
           [ALG.1] [1]: 0.0000089999999772772
           [ALG.2] [1]: 0.0000090000000000000
           [REL.ERR.][ALG.1]: 2.52e-009
           [REL.ERR.][ALG.2]: 0.00e+000

Come previsto, il secondo algoritmo si comporta notevolmente meglio del primo, che compie un errore relativo dell’ordine di circa 10^{−9}.

Esercizio veloce. Testare il codice radicesecgrado.m per l’esempio in cui

p=4.999999999995 \cdot 10^(+4); q=10^(-2); sol=10^(-7);

6.4. Esempio 2: Approssimazione numerica di π. Storicamente sono state scoperte diverse successioni convergenti sempre più rapidamente a π (cf. [18]). In questa sezione ci interesseremo a 3 di queste mostrando sia come le velocità di convergenza possano essere diverse, sia come possano insorgere questioni di instabilità [4, p. 16]. Si implementino le successioni {u_{n}}, {z_{n}}, definite rispettivamente come

\{ \begin{matrix}s_{1} = 1, s_{2} = 1 + \frac{1}{4} \\ u_{1} = 1, u_{2} = 1 + \frac{1}{4} \\ s_{n+1} = s_{n} + \frac{1}{(n+1)^{2}} \\ u_{n+1} = \sqrt{6 s_{n+1}} \end{matrix}

e

\{ \begin{matrix}z_{1} = 1, z_{2} = 2 \\ z_{n+1} = 2^{n− \frac{1}{2}} \sqrt{1 − \sqrt{1 − 4^{1−n} · z_{n}^{2}}} \end{matrix}

che teoricamente convergono a π. Si implementi poi la successione, diciamo {y_{n}}, che si ottiene razionalizzando, cioè moltiplicando numeratore e denominatore per

\sqrt{1 + \sqrt{1 − 4^{1−n} · z_{n}^{2}}}

e si calcolino u_{m}, z_{m} e y_{m} per m = 2, 3, …, 40 (che teoricamente dovrebbero approssimare π). Si disegni in un unico grafico l’andamento dell’errore relativo di u_{n}, z_{n} e y_{n} rispetto a π. A tal proposito ci si aiuti con l’help di Matlab relativo al comando semilogy. I grafici devono avere colori o patterns diversi.

Facoltativo: in un riquadro mettere un legame tra colore (e/o pattern) e successione, usando il comando Matlab legend (aiutarsi con l’help). Ricordiamo che tale comando non esiste in vecchie versioni di GNU Octave. Un esempio del suo utilizzo in Octave è (cf. [7]

legend ("sin (x)");

In seguito scriviamo un’implementazione di quanto richiesto commentando i risultati. Si salvi in un file pigreco.m il codice

% SEQUENZE CONVERGENTI "PI GRECO".
% METODO 1.
s(1)=1; u(1)=1;
s(2)=1.25; u(2)=s(2);
for n=2:40
     s(n+1)=s(n)+(n+1)^(-2);
     u(n+1)=sqrt(6*s(n+1));
     fprintf('\n \t [SEQ.1][INDEX]: %3.0f', n);
     fprintf(' [REL.ERR]: %2.2e', abs(u(n+1)-pi)/pi);
end
rel_err_u=abs(u-pi)/pi;

fprintf('\n');
% METODO 2.
format long
z(1)=1;
z(2)=2;
for n=2:40
     c=(4^(1-n)) * (z(n))^2; inner_sqrt=sqrt(1-c);
     z(n+1)=(2^(n-0.5))*sqrt( 1-inner_sqrt );
     fprintf('\n \t [SEQ.2][N]: %3.0f', n);
     fprintf(' [REL.ERR]: %2.2e', abs(z(n+1)-pi)/pi);
end

Figura 6.1. Grafico che illustra le 3 successioni, rappresentate rispettivamente da ., + e o.

rel_err_z=abs(z-pi)/pi;
fprintf('\n');
% METODO 3.
y(1)=1;
y(2)=2;
for n=2:40
     num=(2^(1/2)) * abs(y(n));
     c=(4^(1-n)) * (z(n))^2;
     inner_sqrt=sqrt(1-c);
     den=sqrt( 1+inner_sqrt );
     y(n+1)=num/den;
     fprintf('\n \t [SEQ.3][N]: %3.0f',n);
     fprintf(' [REL.ERR]: %2.2e', abs(y(n+1)-pi)/pi);
end
rel_err_y=abs(y-pi)/pi;

% SEMILOGY PLOT.
hold on;
semilogy(1:length(u),rel_err_u,'k.');
semilogy(1:length(z),rel_err_z,'m+');
semilogy(1:length(y),rel_err_y,'ro');
hold off;

1. Non tutti i programmi sono functions. Alcuni sono semplicemente un codice che viene interpretato da Matlab (il cosidetto programma principale). Usiamo funzioni solo quando vogliamo introdurre porzioni di codice che possono tornare utili a più

programmi principali, o che semplificano la loro lettura.

2. Piú assegnazioni possono essere scritte in una stessa riga. Per convincersene si osservi la prima riga del file pigreco.m dopo i commenti.

3. L’istruzione descritta dopo for n=2:40 non richiede l’incremento della variabile n, in quanto ciò è automatico.

4. Non serve il ; dopo l’end che chiude il ciclo for.

5. Con un po’ di tecnica il ciclo for è sostituibile col un ciclo while. In questo caso però bisogna incrementare la variabile n.

6. Nel denominatore riguardante l’errore relativo scriviamo pi e non abs(pi) in quanto π = |π|.

7. Nel comando fprintf, utilizziamo \n che manda a capo e \t che esegue un tab (uno spazietto in avanti). Se non si scrive \n, Matlab scriverà di seguito l’output sulla stessa riga, rendendo difficile la lettura dell’errore relativo.

8. Nel comando fprintf alcune parti vengono scritte come testo altre come contenuto di una variabile. Consideriamo ad esempio gli fprintf nella prima successione:

fprintf('\n \t [SEQ.1][INDEX]: \%3.0f',n+1);
fprintf('\n \t [REL.ERR.]: \%2.2e \n', relerru(n+1) );

Dopo essere andati a capo e spaziato a destra, Matlab scrive sul monitor

{\tt {[SEQ.1][INDEX]: } }

e quindi valuta n+1 che viene stampato con tre cifre decimali prima della virgola in notazione decimale. Scrive poi sul monitor [REL.ERR.]:, e accede alla cella di memoria della variabile relerru di cui considera la componente n + 1-sima. Di seguito stampa su monitor il suo valore con due cifre decimali prima della virgola, due cifre decimali dopo la virgola in notazione esponenziale.

9. Il comando semilogy ha come primo argomento l’ascissa (che in questo caso sono gli indici di iterazione 1:41) e quale secondo argomento l’ordinata relerru. Nel grafico (in scala logaritmica nelle ordinate y), vengono disegnate l’i-sima componente dell’ascissa e l’i-sima componente delle ordinate per i = 1, …, dim( relerru). Il grafico viene tenuto grazie al comando di hold on e di seguito si ripete il procedimento per relerrz e relerry. Si osservi che i vettori ascissa e ordinata devono essere della stessa dimensione e dello stesso tipo (cioè entrambi vettori riga o entrambi vettori colonna).

10. Dall’help di semilogy si evince che se la variabile ascissa non viene scritta allora Matlab indicizza automaticamente col vettore di naturali da 1 alla dimensione del vettore ordinata.

Per il risultato del plot si consideri la prima figura. Abbiamo indicato la prima successione con ., la seconda con + e la terza successione con o. Osserviamo che in Octave, i grafici vengono rappresentati differentemente in quanto invece di un cerchietto viene talvolta disegnato un rombo. Dal punto di vista dell’analisi numerica si vede che

1. La prima successione converge molto lentamente a π, la seconda diverge mentre la

terza converge velocemente a π.

2. Per alcuni valori {z_{n}} e {y_{n}} coincidono per alcune iterazioni per poi rispettivamente divergere e convergere a π. Tutto ciò è naturale poichè le due sequenze sono analiticamente (ma non numericamente) equivalenti.

3. Dal grafico dell’errore relativo, la terza successione, dopo aver raggiunto errori relativi prossimi alla precisione di macchina, si assesta ad un errore relativo di circa 10^{−15} (probabilmente per questioni di arrotondamento).

6.7. Esempio 3: Successione ricorrente. In questa sezione mostriamo come alcune formule di ricorrenza in avanti possano essere instabili, mentre d’altro canto le relative versioni all’indietro possono essere stabili [8, Esercizio 1.9, p.35]. Problemi simili con una precisa discussione della propagazione dell’errore sono trattati pure in [4, p. 23, problema 11]

Sia

I_{n} = e^{−1} ∫_{0}^{1} x^{n} e^{x} dx  (6.6)

Per n = 0 si ha

I_{0} = e^{−1} ∫_{0}^{1} e^{x} dx = e^{−1}(e^{1} − 1).

Per n = 1, è facile verificare che, essendo

∫ x exp(x)dx = (x − 1) exp(x) + c,

dal secondo teorema del calcolo integrale (cf. [14]) si ha che I_{1} = e^{−1} e più in generale integrando per parti che

I_{n+1} = e^{−1} (x^{n+1} e^{x} |_{0}^{1} −(n + 1) ∫_{0}^{1} x^{n} e^{x} dx) = 1 − (n + 1) I_{n}.  (6.7)

Da (6.6) essendo l’integranda x^{n} exp x > 0 per x ∈ (0, 1] si ha

I_{n} > 0.

Inoltre x ≤ 1 implica x^{2} ≤ x e più in generale x^{n+1} ≤ x^{n} da cui x^{n+1} exp(x) ≤ x^{n} exp(x) e quindi

I_{n+1} ≤ I_{n}.

La successione I_{n} è positiva e non crescente e quindi ammette limite L finito. Da (6.7), calcolando il limite L per n → ∞ ad ambo i membri si ottiene portando 1 a primo membro

L − 1 = − lim_{n}(n + 1) I_{n}.

L’unica possibilità affinchè lim_{n}(n + 1) I_{n} esista e sia finito è che sia

L = lim_{n} I_{n} = 0.

Infatti se il limite L fosse non nullo, si avrebbe

L − 1 = lim_{n}(n + 1) I_{n} = lim_{n}(n + 1) lim_{n} I_{n} = +∞,

il che è assurdo essendo L < +∞ e quindi L − 1 < +∞.

1. Si calcoli I_{n} per n = 1, …, 99 mediante la successione in avanti

\{ \begin{matrix}s_{1} = e^{−1} \\ s_{n+1} = 1 − (n + 1) s_{n} \end{matrix}

2. Fissato m = 500, si calcoli la successione all’indietro {t_{n}}_{n=1,...,100} definita come

\{ \begin{matrix}t_{2m} = 0 \\ t_{n−1} = \frac{1−t_{n}}{n} \end{matrix}

Si osservi che per raggiungere tale obiettivo bisogna calcolare i termini

t_{m}, t_{m−1}, …, t_{100}, t_{99}, …, t_{2}, t_{1}.

Nota: la successione all’indietro deriva dall’osservare che per n sufficientemente grande I_{n} ≈ 0. Quindi, posto per n sufficientemente grande I_{n} = 0 riscrivendo la successione in avanti come successione all’indietro (cioè I_{n} in funzione di I_{n+1}), abbiamo

I_{n} = \frac{1 − I_{n+1}}{n + 1} .

3. Si disegni in un unico grafico semi-logaritmico (usare semilogy e subplot) l’andamento di |s_{n}| e |t_{n}| per n = 1, …, 100.

4. Si calcoli il valore di t_{1} per diversi valori di m, da m = 1 a m = 10 e lo si confronti con I_{1}.

5. (Esercizio non banale). Si calcolino a mano e_{m}^{(s)} := I_{m} − s_{m} in funzione di e_{1}^{(s)} e e_{1}^{(t)} := I_{1} − t_{1} in funzione di e_{m}^{(t)}. Infine si spieghi l’andamento oscillante della successione {s_{n}}. Suggerimento: si scriva e_{1}^{(s)} = δ (o equivalentemente s_{1} = e^{−1} + δ) ed in seguito si esprima e_{m}^{(s)} := I_{m} − s_{m} in funzione di δ. Da questi conti si vede il motivo del cattivo comportamento della successione in avanti: un piccolo errore sul dato iniziale, ha effetti catastrofici sul risultato al crescere di m.

Di seguito vediamo un’implementazione di quanto richiesto. L’ultimo punto del problema lo lasciamo al lettore. Scriviamo il codice in un file succricorrente.m:

% SUCCESSIONE RICORRENTE.
% SUCCESSIONE "s_n".
s(1)=exp(-1);
for n=1:99

     s(n+1)=1-(n+1)*s(n);
end
% SUCCESSIONE "t_n".
m=500; M=2*m;
t=zeros(M,1); % INIZIALIZZAZIONE "t".
for n=M:-1:2
     j=n-1;
     t(j)=(1-t(n))/n;
end
% PLOT SEMI-LOGARITMICO.
subplot(2,1,1)
semilogy(1:length(s),abs(s),'k-'); hold on;
semilogy(1:length(s),abs(t(1:length(s))),'m-');
hold off;
% ANALISI DI t(1) PER VALORI DIFFERENTI DI "m".
t_1_exact=exp(-1);
for m=1:10
     M(m)=2*m;
     t=zeros(M(m),1); % INIZIALIZZAZIONE "t".
     for n=M(m):-1:2  % SI OSSERVI CHE IL CICLO VA DA "M" A 2.
         t(n-1)=(1-t(n))/n;
     end
     val_t_1(m)=t_1_exact-t(1);
     fprintf('\n \t [M]: %2.0f [VAL.]: %10.15f',M(m),val_t_1(m));
end
subplot(2,1,2)
semilogy(M,abs(val_t_1),'k-');

Alcune osservazioni sul codice

• abbiamo usato un comando del tipo for n=M:-1:2 cosicchè il ciclo for parte da M e arriva a 2, sottraendo di volta in volta 1;

• abbiamo inizializzato con un comando del tipo zeros(M,1) il vettore colonna t, che altrimenti non è in memoria (e quindi le componenti più piccole di t(n-1) sarebbero inaccessibili poichè indefinite);

• abbiamo applicato un grafico semilogaritmico poiché quello usuale non dice granchè (sperimentarlo).

• in figura abbiamo plottato le successioni {|s_{i}|}, {|t_{i}|} e non {s_{i}}, {t_{i}}. Questo non è un problema in quanto

s_{i} → 0 ⇔ |s_{i}| → 0

t_{i} → 0 ⇔ |t_{i}| → 0.

Figura 6.2. Grafico che illustra i valori assoluti assunti dalla successione in avanti (in nero) e all’indietro (in rosa magenta)

D’altro canto la prima successione diverge assumendo anche valori negativi, rendendo difficile la comprensione del grafico.

• Il comando subplot permette di plottare più grafici nella stessa finestra di Matlab/Octave, senza sovrapporli. Per ulteriori informazioni, dall’help di Matlab:

         SUBPLOT Create axes in tiled positions.
            H = SUBPLOT(m,n,p), or SUBPLOT(mnp), breaks the Figure window
            into an m-by-n matrix of small axes, selects the p-th axes for
            for the current plot, and returns the axis handle.  The axes
            are counted along the top row of the Figure window, then the
            second row, etc.  For example,
                 SUBPLOT(2,1,1), PLOT(income)
                 SUBPLOT(2,1,2), PLOT(outgo)
            plots income on the top half of the window and outgo on the
            bottom half.

Numericamente, dopo aver digitato sulla shell di Matlab/Octave succricorrente otteniamo con un po’ di attesa due grafici e i seguenti risultati. Il primo grafico (si veda la Figura 4) mostra un confronto tra i risultati della successione in avanti (quella plottata in nero) e quelli della successione all’indietro (in rosa magenta). Dai ragionamenti qualitativi, è evidente che la prima non fornisce la soluzione, mentre la seconda sembra convergere a 0. Il secondo grafico (si veda la Figura 5) per vari M = 2 · m mostra il comportamento del metodo all’indietro nell’appossimare il valore iniziale exp −1. Evidentemente, al crescere di m, l’approssimazione diventa sempre più accurata. In particolare l’errore assoluto compiuto è stampato nel display dal programma succricorrente.m come segue:

>> succricorrente

Figura 6.3. Grafico che illustra il valore assoluto della successione all’indietro nell’approssimare exp (−1).

           [M]:  2 [VAL.]: 0.132120558828558
           [M]:  4 [VAL.]: 0.007120558828558
           [M]:  6 [VAL.]: 0.000176114384113
           [M]:  8 [VAL.]: 0.000002503273002
           [M]: 10 [VAL.]: 0.000000023114272
           [M]: 12 [VAL.]: 0.000000000149839
           [M]: 14 [VAL.]: 0.000000000000720
           [M]: 16 [VAL.]: 0.000000000000003
           [M]: 18 [VAL.]: 0.000000000000000
           [M]: 20 [VAL.]: 0.000000000000000
>>

7. Complessità computazionale. Ricordato che 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 finale (risultati) in un tempo finito, viene naturale richiedere che esso sia adatto a risolvere una classe di problemi e che sia ottimale rispetto al tempo, al numero di operazioni, alla stabilità e all’accuratezza. Ovviamente ad un maggior numero di operazioni non è detto che corrisponda una miglior accuratezza dei risultati. La complessità computazionale conta il numero di operazioni necessarie per fornire l’approssimazione desiderata.

7.1. Esempio: algoritmo di Horner.. Ci poniamo il problema di valutare il polinomio

p(x) = a_{0} + a_{1} · x + … + a_{n} · x^{n}  (7.1)

in un punto x.

Osserviamo che

p(x) = a_{0} + x · (a_{1} + x · (a_{2} + … + x · (a_{n−1} + x · a_{n})))  (7.2)

Supponiamo sia a = (a_{0}, …, a_{n}) il vettore di dimensione n + 1 delle componenti del polinomio. Possiamo valutare il polinomio tramite i seguenti due algoritmi, il primo che valuta direttamente il polinomio secondo quanto descritto in (7.1), il secondo che effettua la stessa operazione come descritto in (7.2) calcolando dapprima s_{1} = a_{n−1} + x · a_{n}, poi s_{2} = a_{n−2} + x · s_{1} e così via. In Matlab avremo allora

function s=algoritmo1(a,x)
xk=1;
for i=2:length(a)
   xk=xk*x;
   s=s+a(i)*xk;
end

e

L=length(a);
s=a(L); % COMPONENTE a_n IMMAGAZZINATA IN a(n+1).
for i=L-1:-1:1
       s=a(i)+x*s;
end

Se lanciamo il codice demo horner per la valutazione di p(x) = 1+2·x+3·x^{2}+4·x^{3} in x = π

clear all;
a=[1 2 3 4];
x=pi;

y1=algoritmo1(a,x);
y2=algoritmo2(a,x);
format long;
y1
y2

otteniamo

>> demo_horner
ans =
      1.609171052316469e+02
y1 =
      1.609171052316469e+02
y2 =
      1.609171052316469e+02
>>

La differenza sta nella complessità computazionale e non nel risultato numerico. Il primo codice richiede 2n moltiplicazioni e n somme, mentre il secondo algoritmo n moltiplicazioni e n somme.

8. Facoltativo: Un esempio sulla soluzione delle equazioni di secondo grado. Consideriamo l’equazione di secondo grado

x^{2} − 2\sqrt{3}x + 3 = 0  (8.1)

Si nota subito che

x^{2} − 2\sqrt{3}x + 3 = (x − \sqrt{3})^{2}  (8.2)

e quindi (8.2) ha un’unica soluzione \sqrt{3} con molteplicità 2, poichè la derivata di x^{2}−2\sqrt{3}x+3 è la funzione 2x − 2\sqrt{3}x si annulla in \sqrt{3} (cf. [2, p.420]).

Per quanto riguarda il condizionamento di una radice α (cf. [2, p.420]), si prova che se α(ϵ) è la radice di

P_{ϵ}(z) = P (z) + ϵg(z), ϵ > 0  (8.3)

allora se α è semplice

α(ϵ) ≈ α + ϵ (\frac{−g(α)}{P^{1}(α)})  (8.4)

mentre se α è multipla con moltiplicità m, cioè

P (α) = … = P^{(m−1)}(α) = 0 P^{(m)}(α) ≠ 0

allora esiste una radice α(ϵ) di P_{ϵ} tale che

α(ϵ) ≈ α + ϵ^{1/m} (\frac{−m!g(α)}{P^{(m)}(α)})  (8.5)

Approssimiamo in (8.2) il valore \sqrt{3} con 1.7321 fornito da Matlab usando il cosidetto format short. In questo caso, posto g(z) = 2z,

ϵ = −(\sqrt{3} − 1.7321) ≈ 4.9192 · 10^{−5} ≈ 5 · 10^{−5}

abbiamo che

α(ϵ) ≈ α + ϵ^{1/m} (\frac{−m!α}{P^{(2)}(α)})  (8.6)

ed essendo P^{(2)}(z) = 2, m = 2 (molteplicità della radice α = \sqrt{3}), ricaviamo

α(ϵ) ≈ α + (5 · 10^{−5})^{1/2} (\frac{−2!2\sqrt{3}}{2})  (8.7)

in cui

|(5 · 10^{−5})^{1/2} (\frac{−2!2\sqrt{3}}{2})^{1/2} | ≈ 0.02449489742783

Quindi ci si aspetta che una delle radici di

P_{ϵ}(z) = x^{2} − 2 · 1.7321 · x + 3

disti circa 0.02449489742783 da \sqrt{3}.

In effetti, utilizzando semplici modifiche del programma radicisecondogrado, si hanno due radici,

x_{1} ≈ 1.7451541181241765000, x_{2} ≈ 1.7190458818758234000.

ed è , come si vede facilmente dalla shell di Matlab/Octave:

>> 1.7451541181241765000-sqrt(3)
ans =
0.0131
>> 1.7190458818758234000-sqrt(3)
ans =
-0.0130
>>

9. Esercizio. Date le funzioni

f_{1}(x) = 1 − \sqrt{1 − x^{2}}
f_{2}(x) = 1 − x

si calcoli analiticamente il condizionamento

K_{1}(x) = \frac{|x · f′_{1}(x)|}{|f_{1}(x)|}

K_{2}(x) = \frac{|x · f′_{2}(x)|}{|f_{2}(x)|}

Utilizzando Matlab/Octave, si plottino in scala semilogaritmica i grafici di K_{1} e K_{2} nel set di punti x = −1 + 10^{(−3)}, −1 + 2 · 10^{(−3)}, …, 1 − 2 · 10^{(−3)}, 1 − 10^{(−3)},. Si può dire da questi dove la funzione è malcondizionata (cioè assume valori di condizionamento alti)?

10. Alcuni esercizi.

1. Esercizio facile. Si calcoli l’errore relativo tra 1 e il valore che si ottiene valutando in Matlab/Octave

\frac{(1 + η) − 1}{η}

con η = 10^{−1}, 10^{−2}, …, 10^{−15}. Si consideri poi η = 8.8817841970012523E − 16 e si calcoli la medesima quantitá, giustificando i risultati ottenuti [4, p. 5, problema 3]. Suggerimento: quanto vale eps?

2. Esercizio facile. Siano f := tan e g := arctan. Si consideri la funzione composta

f · g(x) := f(g(x)) = tan(arctan(x)) = x, x ∈ (−∞, +∞).

Si calcolino

x = 10^{k}, per k = −20 + 40h, e h = 0, 0.01, 0.02, …, 1

e si valuti l’errore relativo compiuto. Si puó dire che anche numericamente tan(arctan(x)) = x?

Figura 10.1. Grafico che illustra le suddivisioni considerate dall’algoritmo di Archimede, per p = 3, p = 4.

3. Esercizio facoltativo di media difficoltà . Esistono più successioni che approssimano π e sono attribuite ad Archimede. Ad esempio, osservato che questi non è altro che l’area del cerchio avente raggio 1, si inscrivono nel cerchio al variare di p = 2, 3, … dei poligoni regolari aventi 2^{p} lati. Osserviamo che per p = 2, si deve determinare l’area del quadrato inscritto, che con facili conti è uguale a 2. Nel caso generale, si osserva che il poligono consiste di 2^{p} triangoli in cui un lato corrisponde ad un lato del poligono e un vertice con il centro del cerchio (si confronti con le figure). Se indichiamo con θ l’angolo al centro di ogni triangolo, si verifica che la sua area è \frac{sin (θ)}{2}. In particolare per un poligono di 2^{p} lati, si ha θ_{p} = \frac{2π}{2^{p}} e l’area del poligono inscritto, costituito da 2^{p} triangoli uguali, è

I_{p} = 2^{p} \frac{sin (θ_{p})}{2} = 2^{p} \frac{sin (\frac{2π}{2^{p}})}{2}.

In realtà Archimede usò questa ed altre idee per cercare di approssimare π con stime dall’alto e dal basso, via anche poligoni circoscritti. Come citato in [13]: Nel breve lavoro La misura del cerchio viene dimostrato anzitutto che un cerchio equivalente a un triangolo con base eguale alla circonferenza e altezza eguale al raggio. Tale risultato è ottenuto approssimando arbitrariamente il cerchio, dall’interno e dall’esterno, con poligoni regolari inscritti e circoscritti. Con lo stesso procedimento Archimede espone un metodo con il quale può approssimare arbitrariamente il rapporto tra circonferenza e diametro di un cerchio dato, rapporto che oggi si indica con π. Le stime esplicitamente ottenute limitano questo valore fra 22/7 e 3 + 10/71. Secondo [19], Archimede usò un altro algoritmo per approssimare π. Visto che 2π non è altro che la lunghezza della circonferenza del cerchio avente raggio uguale a 1 e questa è maggiore del perimetro a_{k} di un qualsiasi poligono regolare di n = 6 · 2^{k} lati inscritto e minore del perimetro di un qualsiasi poligono regolare di n = 6 · 2^{k} lati circoscritto b_{k}, misurando a_{k} e b_{k} si può affermare che

a_{k} ≤ 2π ≤ b_{k}.

Si può inoltre provare che

a_{k} = 2n sin (\frac{π}{n}),

b_{k} = 2n tan (\frac{π}{n}).

Figura 10.2. Un francobollo raffigurante Archimede (287ac-212ac).

Stimare π come

a_{k} ≤ 2 π ≤ b_{k}.

Per quale n si riesce a determinare π con 10 cifre decimali esatte? Pensare se sia giusto dire b_{k} − a_{k} < 10^{−10} allora π è approssimato da a_{k} con 10 cifre decimali esatte. A. Si approssimi π usando il primo algoritmo di Archimede (basato sull’area di poligoni regolari inscritti) per p = 2 : 20. Si stampi l’errore assoluto e relativo compiuto. B. Stimare π come

a_{k} ≤ π ≤ b_{k}

mediante l’algoritmo di Archimede basato sulla valutazione dei perimetri a_{k}, b_{k} di alcuni poligoni regolari. Per quale n si riesce a determinare π con 10 cifre decimali esatte? Per k = 0, …, 4 si ha

3.00000 ≤ π ≤ 3.46410  (10.1)
3.10583 ≤ π ≤ 3.21539  (10.2)
3.13263 ≤ π ≤ 3.15966  (10.3)
3.13935 ≤ π ≤ 3.14609  (10.4)
3.14103 ≤ π ≤ 3.14271  (10.5)

Osservazione: il fatto che ci sia un π nei secondi membri di alcune espressioni non è un gatto che si morde la coda. Infatti è solo una trascrizione in radianti di un angolo che ha un significato geometrico indipendente dalla quantità π. In questi esercizi, per semplificare la questione si usa però π per calcolare π, cosa che dimostra una volta in più l’abilità di Archimede, che calcolò tali quantità con metodi elementari e geometrici. Si cita il fatto che per k = 4 Archimede riusci’ ad affermare che

\frac{223}{71} ≤ π ≤ \frac{22}{7}.

4. Esercizio, facile facoltativo. Posto h = 10^{−1}, 10^{−2}, …, 10^{−15}, si approssimi la derivata di exp (1) con il rapporto incrementale

\frac{exp (1 + h) − exp (1)}{h}

e quindi si valuti l’errore assoluto compiuto (rispetto alla soluzione exp (1)). L’approssimazione migliora al diminuire di h?

5. Esercizio, facile facoltativo. Si esegua una tabella in cui si studino a due a due le somme, prodotti, divisioni di NaN, Inf, 0, 1 e confrontarle con i noti risultati di indeterminazione della teoria dei limiti.

6. Esercizio, facile facoltativo. Fissato θ = 0 : (2π/1001) : 2π, si valuti sin^{2} (θ) + cos^{2} (θ) e si calcoli l’errore assoluto rispetto il valore 1.

11. Alcune frasi celebri. In [12], oltre a varie curiosità su π si cita un curioso aneddoto riguardante il noto matematico J. H. Conway:

“Un giorno decisi di imparare a memoria le prime mille cifre del pi greco - ricorda Conway - stimolato da mia moglie Larissa, una matematica di origine russa, che aveva bisogno del valore di pi e non ricordava che 3,14. Le insegnai le prime cento cifre che ricordavo già a memoria. Ma questo a lei non bastava e, visto che anch’io non sapevo andare oltre, decidemmo insieme di programmare lo studio di cento nuove cifre ogni giorno, per arrivare almeno a mille, da imparare nei momenti in cui eravamo insieme, al di fuori del nostro lavoro”.

“E’ stato divertente - continua Conway - perchè ogni domenica facevamo una passeggiata fino a Grantchester, una graziosa, piccola cittadina vicino a Cambridge e lungo il percorso recitavamo a turno i gruppi successivi di 20 cifre del pi, come fossero piccole poesie. Venti cifre io e venti cifre mia moglie e così di seguito, alternandoci nella recita: in questo modo siamo arrivati a memorizzare le mille cifre del pi greco”.

Pure la morte di Archimede è coperta dalla leggenda [13]:

Ad un tratto entrò nella stanza un soldato e gli ordinò di andare con lui da Marcello. Archimede rispose che sarebbe andato dopo aver risolto il problema e messa in ordine la dimostrazione. Il soldato si adirò , sguainò la spada e lo uccise. Frasi celebri attribuite ad Archimede

• Noli turbare circulos meos cioè Non turbare i miei cerchi.

• le sue ultime parole pare furono Soldato, stia lontano dal mio disegno.

12. Online. Alcuni siti utili per la comprensione della lezione:

1. http://it.wikipedia.org/wiki/Archimede
2. http://utenti.quipo.it/base5/numeri/pigreco.htm
3. http://it.wikipedia.org/wiki/Pi greco
4. http://it.wikipedia.org/wiki/Teorema fondamentale del calcolo integrale
5. http://it.wikipedia.org/wiki/Floating point
6. http://it.wikipedia.org/wiki/IEEE 754r
7. http://it.wikipedia.org/wiki/IEEE 754
8. http://mathworld.wolfram.com/ArchimedesAlgorithm.html

RIFERIMENTI BIBLIOGRAFICI

[1] K. Atkinson, An Introduction to Numerical Analysis, Wiley, (1989).

[2] V. Comincioli, Analisi Numerica, metodi modelli applicazioni, Mc Graw-Hill, 1990.

[3] V. Comincioli, Problemi di analisi numerica, Mc Graw-Hill, 1991.

[4] M. Frontini e E. Sormani, Fondamenti di Calcolo Numerico, problemi in laboratorio, Apogeo, 2005.

[5] T. Huckle , http://wwwzenger.informatik.tu-muenchen.de/persons/huckle/bugse.html .

[6] J.H. Mathews e K.D. Fink, Numerical Methods using Matlab, Prentice Hall, 1999.

[7] Network Theory, GNU Octave Manual Version, http://www.network-theory.co.uk/docs/octave3/octave 160.html .

[8] A. Quarteroni e F. Saleri, Introduzione al calcolo scientifico, Springer Verlag, 2006.

[9] A. Quarteroni, R. Sacco e F. Saleri, Matematica Numerica, Springer Verlag, 1998.

[10] J. Stoer, Introduzione all’analisi numerica, Zanichelli, 1984.

[11] The MathWorks Inc., Numerical Computing with Matlab,http://www.mathworks.com/moler .

[12] Web Page , http://utenti.quipo.it/base5/numeri/pigreco.htm .

[13] Wikipedia (Archimede), http://it.wikipedia.org/wiki/Archimede

[14] Wikipedia (Teorema Fondamentale del Calcolo Integrale), http://it.wikipedia.org/wiki/Teorema fondamentale del calcolo integrale

[15] Wikipedia (Floating Point), http://it.wikipedia.org/wiki/Floating point

[16] Wikipedia (IEEE 754), http://it.wikipedia.org/wiki/IEEE 754

[17] Wikipedia (IEEE 754r), http://it.wikipedia.org/wiki/IEEE 754r

[18] Wikipedia (Pi Greco), http://it.wikipedia.org/wiki/Pi greco

[19] Wolfram (Archimedes’ Algorithm), http://mathworld.wolfram.com/ArchimedesAlgorithm.html