Sulla propagazione degli errori

3 maggio 2008

1 Esercitazione sulla propagazione degli errori

1.1 Facoltativo: Rappresentazione dei numeri al calcolatore

Risolvere un problema mediante calcolo numerico ed al calcolatore (il calcolo numerico era nato come calcolo a mano, ma oggi questo modus operandi è praticamente caduto in disuso, dato il bassissimo costo del calcolo automatico), significa utilizzare un’aritmetica a precisione finita, e quindi produrre soluzioni approssimate e non esatte. Soluzioni esatte possono essere prodotte se si riesce a risolvere il problema mediante calcolo simbolico.

Per quanto concerne la rappresentazione numerica, i calcolatori moderni rispettano in genere lo standard IEEE 754 che andiamo brevemente a descrivere.

Si intende per bit una cifra nel sistema (usualmente binario) usato [1, p. 48], [7, p. 47].

Nello standard IEEE 754, i numeri in virgola mobile a precisione singola occupano una parola da 32 bits, mentre quelli a precisione doppia occupano due parole consecutive da 32 bits. Un numero non nullo-normalizzato si scrive come

( 1)^{}2^{}(1 )  (1)

dove 0 oppure 1 determina il segno, è l’esponente cui va sottratto il valore del BIAS ed la mantissa. La scrittura 1 significa 1 dove _{1}^{}_{}2^{}, con _{}0 oppure _{}1. Alcuni esempi

Precisione Ebits Bias FBits
Singola 8 127 23
Doppia  11 1023 52

Per capire meglio tale notazione consideriamo il link

          http://babbage.cs.qc.edu/IEEE-754/Decimal.html

in cui viene convertito un numero decimale in forma binaria, secondo lo standard IEEE 754. Partiamo con un esempio semplice. Il numero 12 corrisponde in singola precisione a

[0 10000010 10000000000000000000000]

Cerchiamo di comprenderne il significato:

Propagazione degli errori

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!) 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^{}2^{3} 8.

3. 10000000000000000000000 è la componente relativa alla mantissa. La quantità , letta da sinistra a destra, riguarda la parte decimale in (1). Quindi

1 2^{1} 0 2^{2}  0 5. Ora (1 ) si legge come 1  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!) 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^{}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. 11111000000000000000000 è la componente relativa alla mantissa. La quantità , letta da sinistra a destra, riguarda la parte decimale in (1). Quindi 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 ) si legge come 1 1 96875.

In altre parole la rappresentazione di 126 è esatta ed il numero è rappresentato come 126 64 1 96875.

Dalla tabella riguardante la precisione singola e i bit, 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. Infatti, la distanza tra un numero in virgola mobile ed il successivo, è :

2-nbitsF 2E-bias

La sua cardinalità è piuttosto elevata: in doppia precisione i numeri macchina sono

2^{6}4 1 844674407370955 019

Si dice precisione di macchina, quel numero che rappresenta la distanza tra 1 ed il successivo numero in virgola mobile. In Matlab tale valore è rappresentato da eps. Vediamo di calcolarlo. Si scriva sulla shell di Linux

                      pico myeps.m

e nel file il codice Matlab

myeps=1;
while 1+myeps > 1;
      myeps=myeps/2;
end
myeps=2*myeps;
format long e
myeps-eps

Il codice Matlab parte con myeps=1 e continua a dividere per 2 finché trova un numero tale che sommato a 1 dà (numericamente, non analiticamente) 1. Vista la definizione di eps alla fine bisogna porre myeps 2 myeps (altrimenti myeps eps). E infatti myeps-eps 0.

1.2 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 vari-

abile realmin. Il limite superiore è pari a 2^{2Ebits 2} (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 225073858507201 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, Ebits 11, ed Fbits 52. Conseguentemente

realmax  2^{2Ebits 2} (1 (1 2-Fbits))
2^{2^{11} 2 1023} (1 (1 2^{52}))
2^{2^{11} 2 1023} (2 2^{52})
1 797693134862316 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*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^{5}2) 4 940656458412465e 324

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

( 1)^{}2^{}(0 )  (2)

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*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 .* (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.
>>

Lo standard richiede che il risultato delle operazioni di addizione, sottrazione, moltiplicazione e divisione sia arrotondato esattamente, cioè il risultato deve essere calcolato esattamente e poi arrotondato al numero in virgola mobile più vicino.

Problemi in tal senso provengono dalla sottrazione. Per questo motivo, nei microprocessori moderni la sottrazione viene effettuata utilizzando la tecnica bit di guardia. Non tratteremo questo punto dettagliatamente. Per curiosità si consulti il sito

        http://babbage.cs.qc.edu/IEEE-754/References.xhtml

1.3 Propagazione degli errori di arrotondamento

Ogni calcolatore può rappresentare esattamente solo un numero finito di numeri reali, i cosidetti numeri macchina. Per i rimanenti, seguendo uno standard quale IEEE 754, può fornirne solo un’approssimazione. Di conseguenza, le operazioni aritmetiche di base possono in generale essere soggette ad un errore nel risultato ed è dunque necessario conoscere l’entità di quest’ultimo e l’effetto che può avere in un algoritmo. Questo problema non è da trascurarsi. Si possono trovare in internet vari siti in cui queste approssimazioni, dovute ad esempio al cambiamento di unità di misura, hanno portato a disastri quali la perdita di satelliti o la distruzione di stazioni petrolifere come indicato ad esempio in [4]. Passiamo quindi in dettaglio all’analisi del problema. Sia 0 un numero reale del tipo

( 1)^{}2^{( bias)} (1 )

in cui supponiamo 2^{( nbitsF)}/2. In precisione singola, secondo IEEE 754, Ebits 8, bias 127, Fbits 23. La scrittura 1 significa 1 dove _{1}^{}_{}2^{}, con 0 oppure _{}1.

Allora, la sua rappresentazione in virgola mobile, sarà

fl( ) ( 1)^{}2^{( bias)} 1

con un errore relativo (di arrotondamento)

\frac{fl( )}{}\frac{2^{( bias)}}{2^{( bias)} (1 )} \frac{}{1}

Ricordiamo che fissato un vettore da approssimare

( _{1}^{}_{}^{}) �^{}

si definisce 1. errore assoluto tra e ^{}la quantità

dove

\overline{_{1}^{}_{}^{2}}  ( _{1}  _{});

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

\frac{^{}}{^{}}

Si noti che se ^{}0, cioè ^{}0 per la proprietà delle norme [?, p.481], allora non ha senso la scrittura propria dell’errore relativo. 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( ) è il numero macchina che corrisponde a , denotati con

fl(fl( ) fl( ))  (3)
fl(fl( ) fl( ))  (4)
fl(fl( ) fl( ))  (5)
fl(fl( ) : fl( ))  (6)

e per una delle operazioni precedentemente introdotte, cioè una tra , , , corrispondenti a op cioè una tra , , , :, posto

\frac{fl( )}{}(7)
_{}^{}\frac{( op ) ( )}{}(8)

si ha dopo qualche non banale conto (cf. [1, p.78])

_{}^{}\frac{}{}_{}\frac{}{}_{}(9)
_{}^{}\frac{}{}_{}\frac{}{}_{}(10)
_{}^{}_{}_{}(11)
_{}^{}_{}_{}(12)

il che mostra il pericolo di perdita di accuratezza nella somma e nella sottrazione qualora rispettivamente 0 e 0. Vediamo di seguito alcuni esempi in cui è fondamentale un adeguato trattamento dell’errore per giungere ad una buona approssimazione del risultato esatto.

1.4 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 ^{2} 2 , con \overline{^{2}} 0 calcoliamo la radice

\overline{^{2}}  (13)

Osserviamo che essendo \overline{^{2}} 0 le radici

\overline{^{2}}  (14)

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

\overline{^{2}} \frac{(  \overline{^{2}} )(  \overline{^{2}} )}{(  \overline{^{2}} )} \frac{}{(  \overline{^{2}} )}  (15)

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 il suo comportamento rispetto alla propagazione degli errori è più usuale il termine di stabilità [1, p. 66].

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

1. il primo algoritmo non è numericamente stabile qualora ;

2. il secondo algoritmo è numericamente stabile qualora .

In [8, p.22], si suggerisce un test interessante per

1000 0 018000000081

la cui soluzione esatta è 0 9 10^{5} . Secondo [2], p. 11 è notevole l’esperimento in cui

4 99999999999510^{4} 10^{2}

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

1.4.1 Implementazione in Matlab

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

% p=4.999999999995*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. Testare il codice radicesecgrado.m per l’esempio in cui

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

1.5 Approssimazione numerica di

Storicamente sono state scoperte diverse successioni convergenti sempre più rapidamente a (cf. [16]). 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à [3, p. 16]. 1.5.1 Successioni convergenti a

: tests numerici Si implementino le successioni _{}, _{}, definite rispettivamente come

 ;  ;  ;  ;  ;  ;  ;  ; \begin{matrix}_{1} 1 _{2} 1 \frac{1}{4} \\ _{1} 1 _{2} 1 \frac{1}{4} \\ _{1}  _{}\frac{1}{( 1)^{2}} \\ _{1}  \overline{6 _{1}} \end{matrix}

e

 ;  ; \begin{matrix}_{1} 1 _{2} 2 \\ _{1} 2^{\frac{1}{2}}  \overline{1 \overline{1 4^{1} _{}^{2}}} \end{matrix}

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

\overline{1 \overline{1 4^{1} _{}^{2}}}

e si calcolino _{}, _{}e _{}per 2 3 40. Si disegni in un unico grafico l’andamento dell’errore relativo di _{}, _{}e _{}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 GNU Octave.

1.5.2 Risoluzione

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
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.
semilogy(1:length(u),rel_err_u,'k.'); hold on;
semilogy(1:length(z),rel_err_z,'m+'); hold on;
semilogy(1:length(y),rel_err_y,'ro');

Figure 1: Grafico che illustra le 3 successioni, rappresentate rispettivamente da ., + e o.

1.5.3 Commenti alla Risoluzione

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 , 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 .

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

                      [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 1sima. 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 ), vengono disegnate l’ -sima componente dell’ascissa e l’ -sima componente delle ordinate per 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 fi 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 _{}e _{}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).

Figure 2: Grafico che illustra le suddivisioni considerate dall’algoritmo di Archimede, per 2, 3.

1.6 Esercizio: alcune successioni di Archimede

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 2 3 dei poligoni regolari aventi 2^{}lati. Osserviamo che per 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^{}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^{}lati, si ha _{}\frac{2}{2^{}}e l’area del poligono inscritto, costituito da 2^{}triangoli uguali, è

2^{}\frac{sin( _{})}{2} 2^{}\frac{sin(\frac{2}{2^{}})}{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 [11]:

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 [17], 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 _{}di un qualsiasi poligono regolare di 6 2^{} lati inscritto e minore del perimetro di un qualsiasi poligono regolare di 6 2^{}lati

Figure 3: Un francobollo raffigurante Archimede (287ac-212ac).

circoscritto _{}, misurando _{}e _{}si può affermare che

2 _{}

Si può inoltre provare che 2 tan \frac{}{}

2 sin \frac{}{}

Stimare come

Per quale si riesce a determinare con 10 cifre decimali esatte? Esercizio 1. Si approssimi usando il primo algoritmo di Archimede (basato sull’area di poligoni regolari inscritti) per 2 : 20. Si stampi l’errore assoluto e relativo compiuto. Esercizio 2. Stimare come

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

3 00000  3 46410  (16)
3 10583  3 21539  (17)
3 13263  3 15966  (18)
3 13935  3 14609  (19)
3 14103  3 14271  (20)

Si cita il fatto che per 4 Archimede riusci’ ad affermare che

\frac{223}{71} \frac{22}{7}

1.7 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 [6, Esercizio 1.9, p.35]. Problemi simili con una precisa discussione della propagazione dell’errore sono trattati pure in [3, p. 23, problema 11]

Sia

1 _{0}^{1} ^{}^{}(21)

Per 0 si ha

0 ^{1} _{0}^{1} ^{}^{1} ( ^{1} 1)

Per 1, è facile verificare che, essendo

exp( ) ( 1) exp( )

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

1  ^{1} ^{1} ^{}_{0}^{1} ( 1) _{0}^{1} ^{}^{}1 ( 1) _{}(22)

Da (21) essendo l’integranda ^{}exp 0 per (0 1] si ha

0

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

1 _{}

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

1 lim_{}( 1) _{}

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

lim_{}_{}0

Infatti se il limite fosse non nullo, si avrebbe

1 lim_{}( 1) _{}lim_{}( 1) lim_{}_{}

il che è assurdo essendo e quindi 1 .

1. Si calcoli _{}per 1 99 mediante la successione in avanti

\begin{matrix}_{1} ^{1} \\ _{1} 1 ( 1) _{} \end{matrix}

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

\begin{matrix}_{2} 0 \\ _{1} \frac{1 _{}}{} \end{matrix}

Si osservi che per raggiungere tale obiettivo bisogna calcolare i termini

1  _{100} _{99}  _{2} _{1}

3. Si disegni in un unico grafico semi-logaritmico (usare semilogy e subplot) l’andamento di _{}e _{}per 1 100.

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

5. Si calcolino infine a mano _{}^{( )}: _{}_{}in funzione di _{1}^{( )} e _{1}^{( )}: _{1} _{1} in funzione di _{}^{( )}. Infine si spieghi l’andamento oscillante della successione .

1.7.1 Risoluzione

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);

Figure 4: Grafico che illustra i valori assoluti assunti dalla successione in avanti (in nero) e all’indietro (in rosa magenta)

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-');

1.7.2 Commenti alla Risoluzione

Alcune osservazioni sul codice

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

• abbiamo inizializzato con un comando del tipo zeros(M,1) il vettore colonna , 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 _{}, _{}e non _{}, _{}. Questo non è un problema in quanto

0 _{}0

0 _{}0

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 2) 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 1.2) per vari 2 mostra il comportamento del metodo all’indietro nell’appossimare il valore iniziale exp 1. Evidentemente, al crescere di , l’approssimazione diventa sempre più precisa. In particolare l’errore assoluto compiuto è stampato nel display dal programma succricorrente.m come segue:

>> succricorrente
          [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

Figure 5: Grafico che illustra il valore assoluto della successione all’indietro nell’approssimare exp( 1).

          [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
>>

2 Alcuni esercizi

2.1 Esercizio 1

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 8817841970012523 16 e si calcoli la medesima quantitá, giustificando i risultati ottenuti [3, p. 5, problema 3]. Suggerimento: quanto vale eps?

2.2 Esercizio 2

Siano : tane : arctan. Si consideri la funzione composta

( ) : ( ( )) tan(arctan( ))  (  )

Si calcolino 10^{}per 20 40 e 0 0 01 0 02 1

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

3 Alcune frasi celebri

In [10], 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 [11]:

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.

4 Online

Alcuni siti utili per la comprensione della lezione:

http://it.wikipedia.org/wiki/Archimede
http://utenti.quipo.it/base5/numeri/pigreco.htm
http://it.wikipedia.org/wiki/Pi_greco
http://it.wikipedia.org/wiki/Teorema_fondamentale_del_calcolo_integrale
http://it.wikipedia.org/wiki/Floating_point
http://it.wikipedia.org/wiki/IEEE_754r
http://it.wikipedia.org/wiki/IEEE_754
http://mathworld.wolfram.com/ArchimedesAlgorithm.html

References

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

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

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

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

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

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

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

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

[9] The MathWorks Inc., Numerical Computing with Mat-

   lab,http://www.mathworks.com/moler.
[10] Web Page , http://utenti.quipo.it/base5/numeri/pigreco.htm.
[11] Wikipedia, http://it.wikipedia.org/wiki/Archimede
[12] Wikipedia, http://it.wikipedia.org/wiki/Teorema_fondamentale_del_calco
   lo_integrale
[13] Wikipedia, http://it.wikipedia.org/wiki/Floating_point
[14] Wikipedia, http://it.wikipedia.org/wiki/IEEE_754
[15] Wikipedia, http://it.wikipedia.org/wiki/IEEE_754r
[16] Wikipedia, http://it.wikipedia.org/wiki/Pi_greco
[17] Wikipedia, http://mathworld.wolfram.com/ArchimedesAlgorithm.html