
% Copyright (C) 2007 Alvise Sommariva. This is a free document; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% at your option) any later version. This document is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details. You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA}

\documentclass[final]{siamltex}


%------------------------------------------------------------------------------------------------------------------------
% USEPACKAGE.
%------------------------------------------------------------------------------------------------------------------------

\usepackage[italian]{babel}
\usepackage[dvips]{graphicx}
\usepackage{amssymb}
\usepackage{amsmath}
\usepackage{times}
\usepackage{epsf}
\usepackage{graphicx,psfrag,color,pst-grad}
\usepackage{latexsym}
\usepackage{psboxit,pstricks}
\usepackage{fancyhdr}
\usepackage{fancyvrb}
\usepackage[colorlinks=true,linkcolor=purple,urlcolor=purple]{hyperref}

\usepackage{float}
\floatstyle{ruled}
\newfloat{program}{thp}{lop}
\floatname{program}{Codice}
\pdfpagewidth 8.5in
\pdfpageheight 11in

% \usepackage{palatino}
% \usepackage[latin1]{inputenc}
\usepackage{calc}
\usepackage{multicol}



%------------------------------------------------------------------------------------------------------------------------
% NEW COMMAND.
%------------------------------------------------------------------------------------------------------------------------

\newcommand{\ptenm}[2]{#1\mbox{-E#2}}
\newcommand{\ptenp}[2]{#1\mbox{+E#2}}
\newcommand{\mye}{\`e }
\newcommand{\myer}{\'e }
\newcommand{\myu}{\`u }
\newcommand{\myur}{\'u }
\newcommand{\mya}{\`a }
\newcommand{\myo}{\`o }
\newcommand{\myod}{\'o }
\newcommand{\myi}{\`i  }
\newcommand{\myid}{\'i }

\newcommand{\myes}{\`e}
\newcommand{\myers}{\'e}
\newcommand{\myus}{\`u}
\newcommand{\myurs}{\'u}
\newcommand{\myas}{\`a}
\newcommand{\myos}{\`o}
\newcommand{\myods}{\'o}
\newcommand{\myis}{\`i}
\newcommand{\myids}{\'i}

\newcommand{\myfl}{\mbox{fl}}
\newcommand{\sgn}{\mbox{sign}}
\newcommand{\hsp}{\hspace{0.2cm}}
\newcommand{\smsp}{\vspace{0.25cm}}
\newcommand{\tnsp}{\vspace{0.1cm}}
\newcommand{\Rset}{\mathbb{R}}
\newcommand{\bs}{\boldsymbol}
\newcommand{\mylbar}{|}
\newcommand{\myrbar}{|}

%\pagestyle{fancy}
%\headheight 35pt
%\rhead{}
%\chead{}
%\lhead{}
%\rfoot{}
%\cfoot{\thepage}
%\lfoot{}



%------------------------------------------------------------------------------------------------------------------------
% NEW THEOREM.
%------------------------------------------------------------------------------------------------------------------------

%\newtheorem{theorem}{Theorem}[section]
%\newtheorem{lemma}[theorem]{Lemma}
%\newtheorem{proposition}[theorem]{Proposition}
%\newtheorem{corollary}[theorem]{Corollary}
% \newtheorem{remark}[theorem]{Nota}
%\newenvironment{proof}[1][Proof]{\begin{trivlist}
%\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}
%\newenvironment{definition}[1][Definition:]{\begin{trivlist}
%\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}
%\newenvironment{example}[1][Esempio]{\begin{trivlist}
%\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}
%\newenvironment{remark}[1][Nota]{\begin{trivlist}
%\item[\hskip \labelsep {\bfseries #1}]}{\end{trivlist}}


%------------------------------------------------------------------------------------------------------------------------
% NEW COMMAND.
%------------------------------------------------------------------------------------------------------------------------

\renewcommand{\figurename}{Figura}

%\newcommand{\qed}{\nobreak \ifvmode \relax \else
%     \ifdim\lastskip<1.5em \hskip-\lastskip
%    \hskip1.5em plus0em minus0.5em \fi \nobreak
% \vrule height0.75em width0.5em depth0.25em\fi}



%------------------------------------------------------------------------------------------------------------------------
% DEF.
%------------------------------------------------------------------------------------------------------------------------

\def\P{{\mathbb P}}
\def\n{{\mathcal N}}
\def\R{\mathbb{R}}
\def\N{\mathbb{N}}
\def\Z{\mathbb{Z}}
\def\F{{\mathcal F}}
\newcommand{\LL}{\mathcal L}
\def\C{\mathbb{C}}
\def\K{\mathbb{K}}
\def\T{\mathbb{T}}
\def\qed{\vbox{\hrule height0.6pt\hbox{%
   \vrule height1.3ex width0.6pt\hskip0.8ex
   \vrule width0.6pt}\hrule height0.6pt
  }}

\begin{document}

\title{\bf Numeri macchina \thanks{Ultima revisione: \today}}

\author{A. Sommariva\thanks{Dipartimento di Matematica,
Universit\mya degli Studi di  Padova, stanza 419, via Trieste 63, 35121 Padova, Italia ({\tt 
alvise@math.unipd.it}).} Telefono: +39-049-8271350.
}

\maketitle

\begin{knowledge}
Equazioni di secondo grado, successioni, formula di Taylor in una e due variabili, integrali di Riemann.
\end{knowledge}

\begin{achieved}
In questa sezione analizziamo la rappresentazione floating point di un numero. Di seguito introduciamo il concetto di errore assoluto e relativo, condizionamento, stabilit\mya e complessit\mya computazionale con esempi.
\end{achieved}

\section{Numeri macchina e loro propriet\myas}
Fissato un numero naturale $\beta > 1$ \mye possibile vedere che ogni numero reale $x \neq 0$ ha una unica rappresentazione del tipo
\begin{equation}\label{E1}
x=\sgn(x)  \beta^e  \sum_{k=1}^{+\infty} d_k \beta^{-k}, \hspace{0.3cm} 0 < d_1 , \, 0\leq d_k \leq \beta-1, \, e\in{\mathbb{Z}}
\end{equation}
dove 
$$\sgn(y)=\left\{ \begin{array}{l}1, \hspace{0.575cm} y > 0 \\ -1, \hspace{0.3cm} y <0 \\ 0, \hspace{0.575cm} y = 0.  \end{array} \right.$$
\mye la funzione segno.

Osserviamo che storicamente, la scelta della base $\beta$ \mye risultata molto eterogenea:
\begin{itemize}
\item la base utilizzata dalla maggior parte dei computer attuali \mye $\beta=2$, 
\item il computer russo SETUN ($\approx 1957$) usava la base $\beta=3$,
\item il PDP-10 ($\approx 1965$), Burroughs 570 e 6700 utilizzavano la base ottale, cio\mye $\beta=8$,
\item la base che usiamo nella vita di tutti i giorni \mye $\beta=10$, ed \mye stata una scelta comune anche per computer del passato come l'HP-11C, HP-15C ($\approx 1980$) o Eniac ($\approx 1947$),
\item esistono popoli in Nigeria, India, Nepal, la cui \mye base $\beta=12$,
\item molti computers, come l'IBM 360, l'IBM 3033  ($\approx 1977$), usavano la base $\beta=16$,
\item i Maya usavano la base $\beta=20$. 
\item in tavolette babilonesi si utilizzava la base $\beta=60$.
\end{itemize}

Dall'avvio dello standard IEEE 754-1985, in cui venne realizzato lo standard per l'aritmetica binaria floating point, la base pi\myu utilizzata risulta $\beta=2$.  \cite{BOth}.

Vediamo quale esempio come rappresentare il numero $\pi$ al variare di due classiche scelte di $\beta$:
\begin{itemize}
\item base 10: $3,14159 26535 89793 23846 26433 83279 50288 41971 69399 \ldots$ secondo la rappresentazione sopra descritta e'
$$
+1 \cdot 10^1 \cdot (3 \cdot 10^{-1} + 1  \cdot 10^{-2} + 4 \cdot 10^{-2} + \ldots)
$$
che viene spesso denotato come
$$
\pi=(+1) \cdot (0.314159 26535 89793 23846 26433 83279 50288 41971 69399 \ldots)_{10} \cdot 10^1;
$$
\item base 2: $11,0010010000111111011010101000100010000101 \ldots$  secondo la rappresentazione sopra descritta e'
$$
+1  \cdot 2^2 \cdot (1 \cdot 2^{-1} + 1  \cdot 2^{-2} + 0 \cdot 2^{-3} \ldots + 0 \cdot 2^{-4} + 1 \cdot 2^{-5} \ldots )=+4(1/2+1/4+1/32+\ldots),
$$
che viene spesso denotato come
$$
\pi=(+1) \cdot (0.110010010000111111011010101000100010000101 \ldots)_2 \cdot 2^2.
$$
\end{itemize}



La particolarit\mya della sommatoria in (\ref{E1}), ovvero
\begin{equation}\label{E1R}
x=\sgn(x)  \beta^e  \sum_{k=1}^{+\infty} d_k \beta^{-k}, \hspace{0.3cm} 0 < d_1 , \, 0\leq d_k \leq \beta-1, \, e\in{\mathbb{Z}}
\end{equation}
 \mye ovviamente tale che $x$ potrebbe essere non rappresentabile esattamente dal calcolatore se esistono infiniti $d_k >0$, in quanto il calcolatore \mye capace di immagazzinare solo un numero finito di tali $d_k$. 
 
 Ad esempio, non \mye difficile vedere che:
\begin{itemize}
\item i numeri irrazionali come $\pi$, indipendentemente dalla base $\beta$, hanno sempre un numero infinito di cifre $d_k$;
\item in base $10$ alcuni numeri razionali hanno comunque un numero infinito di cifre $d_k$ (si pensi ai numeri periodici come $1/3=0.333\ldots$, ma comunque in certe basi possono avere una rappresentazione finita (nel nostro esempio $1/3=1 \cdot 3^{-1}$ per $\beta=3$).
\end{itemize}

\vspace{0.25cm}

In virt\'u di di considerazioni sulla finitezza delle cifre immagazzibili, ipotesi necessaria perch\'e il numero sia utilizzabile da un computer, risulta naturale la seguente

\vspace{0.3cm}

 \begin{definition}
L'insieme dei numeri macchina $F(\beta,t,L,U)$ di $\mathbb{R}$ \mye costituito da quei numeri $y$ che sono $0$ oppure sono rappresentati con il {\sl{sistema a virgola mobile normalizzata o floating point normalizzato}} da $$y=\sgn(y) \cdot (0.d_1 \ldots d_t)_{\beta} \beta^e= \sgn(y) \cdot \beta^e \sum_{k=1}^{t} d_k \beta^{-k}, \hspace{0.3cm} d_1 >0, \, 0\leq d_k \leq \beta-1$$
con 

\vspace{0.3cm}

\begin{itemize}
\item $\beta > 1$ numero naturale detto {\sl{base}}, $t$ numero prefissato di cifre di {\sl{mantissa}},  
\item $e \in \mathbb{Z}$ l'esponente intero tale che $L \leq e \leq U$. 
\end{itemize}

\vspace{0.3cm}

La parola {\sl{normalizzato}} sottolinea che $d_1 > 0$. Si pu\myo provare che $\sum_{k=1}^{t} d_k \beta^{-k}<1$ e quindi $e={\tt{floor}}(1+{\log_{\beta}}(x))$.

Ogni elemento di $F(\beta,t,L,U)$ \mye detto {\sl{numero macchina}}, ed evidentemente \mye un numero razionale. 

\end{definition}

\vspace{0.3cm}

\begin{remark}
Osserviamo che se $e > 0$ allora 
\begin{eqnarray}
y&=&\sgn(y) \cdot (0.d_1 \ldots d_t)_{\beta} \beta^e= \sgn(y) \cdot \beta^e \sum_{k=1}^{t} d_k \beta^{-k} \nonumber \\
&=&\sgn(y) \cdot ( \underbrace{\sum_{k=1}^{e} d_k \beta^{e-k}}_{\mbox{parte intera}} +\underbrace{ \sum_{k=e+1}^{t} d_k \beta^{e-k})}_{\mbox{parte frazionaria}} \nonumber
\end{eqnarray}
dove l'ultima somma \mye nulla per definizione se $t < e+1$.
\end{remark}

\vspace{0.3cm}


%Quale ulteriore esempio vogliamo vedere se $100$ appartiene a $F(10,6,-6,8)$. Essendo l'esponente corretto $e={\tt{floor}}(1+{\log_{10}}(100))=3$, si vede subito che $100=10^3 \cdot 0.1=10^3 \cdot (0.100000) \beta^e$.

\vspace{0.3cm}%

\begin{remark}
Si osservi che il sistema floating-point \mye solo una delle scelte possibili, in quanto una volta si utilizzava un sistema a {\sl{virgola fissa}}, ovvero con un numero fissato di cifre prima della virgola e un certo numero di cifre dopo la virgola, abbandonato (nonostante presentasse alcuni vantaggi) perch\mye l'intervallo di valori rappresentati \mye modesto e la precisione dei numeri frazionari scarsa.
\end{remark}

\vspace{0.3cm}

Di ogni numero macchina si immagazzinano 

\vspace{0.3cm}

\begin{itemize}
\item una cifra per il segno (ovvero $0$ per positivo e $1$ per negativo), 
\item le cifre della mantissa,
\item le cifre dell'esponente. 
\end{itemize}

\vspace{0.3cm}

Ognuna di queste cifre \mye detta {\sl{bit}}. Ogni numero \mye registrato in un vettore detto {\sl{parola}} di un certo numero di $M$ {\sl{bits}}. 

\subsection{Precisione singola e doppia} Descriviamo di seguito i pi\myu comuni insiemi di numeri floating-point. 
Usualmente se $M$ \mye il numero di bits necessari per un numero in precisione {\sl{singola}}, nel caso di precisione {\sl{doppia}} questo \mye usualmente composto da 2 {\sl{parole}} di $M$ bits. 

Classici esempi sono per quanto riguarda i processori che seguono lo standard IEEE754

\vspace{0.3cm}

\begin{itemize}
\item per la precisione singola $F(2,24,-126,127)$ in cui ogni numero macchina occupa 32 bit (o 4 byte essendo un byte pari a 8 bit) 
\item per la doppia $F(2,53,-1022,1023)$ in cui ogni numero macchina occupa $64$ bit. 
\end{itemize}

\vspace{0.3cm}

\begin{remark}
In altri sistemi il numero di cifre per esponente e mantissa sono risultati molteplici, ad esempio su
\begin{itemize}
\item Zuse Z1 ($\approx 1936$), la mantissa corrispondeva a $22$ cifre (14 per la mantissa e 8 per l'esponente),
\item IBM 3033 in precisione doppia si aveva, oltre alla base $\beta=16$, una mantissa di $14$ cifre,
\item PRIME 850, che lavorava in base $\beta=2$, si aveva in precisione singola $23$ cifre di mantissa, mentre in precisione doppia $47$ cifre di mantissa, \cite[p.13]{ATK}.
\end{itemize}
\end{remark}

\vspace{0.25cm}

Vediamo cosa serva immagazzinare in precisione doppia. Essendo un numero non nullo di segno $+1$ o $-1$, basta un solo bit per il segno. Delle 53 cifre della mantissa, solo 52 devono essere ricordate (la prima per forza \mye $1$ in virt\myu della normalizzazione). Infine, per l'esponente, vista una trattazione particolare che tiene conto anche di possibili problemi, servono 11 bit.

Osserviamo che 

\begin{itemize}
\item per $e \geq t$ con $t$ numero di cifre della mantissa, essendo $\beta > 1$ un numero naturale, allora 
$$\pm \beta^e \sum_{k=1}^{t} d_k \beta^{-k}=\pm \sum_{k=1}^{t} d_k \beta^{e-k}$$
\mye un numero intero, essendo $e-k>0$. Di conseguenza i numeri macchina pi\myu grandi in valore assoluto sono interi;

%Si vede facilmente che tutti gli interi in valore assoluto inferiori o uguali a $\beta^U \cdot (1-\beta^{-t})$ sono in particolare nell'insieme $F(\beta,t,L,U)$;

\item l'insieme dei numeri macchina $F(\beta,t,L,U)$ ha un numero finito di elementi e si pu\myo dimostrare che la sua cardinalit\mya \mye
$$2 (\beta-1) \beta^{t-1}(U-L+1)+1.$$
Di conseguenza

\vspace{0.25cm}

\begin{itemize}
\item  in singola precisione $F(2,24,-126,127)$ i numeri sono circa $4.3 \cdot 10^9$;
\item in precisione doppia $F(2,53,-1022,1023)$ sono circa $1.8 \cdot 10^{19}$.
\end{itemize}



\item ponendo $d_1=1$, $d_2=\ldots=d_{t}=0$ ed $e=L$ otteniamo che il pi\myu piccolo numero macchina (e quindi normalizzato!) positivo rappresentabile \mye $\beta^{L-1}$.
\item ponendo $d_k=\beta-1$, per $k=1,\ldots,t$ ed $e=U$ otteniamo che il massimo numero macchina rappresentabile \mye $$\beta^U (1-\beta^{-t});$$
\item essendo un numero macchina del tipo
$$x = \pm \beta^e  \sum_{k=1}^{t} d_k \beta^{-k} $$
con $d_1 >0$ e $0 \leq d_k \leq \beta-1$, da
$$
\sum_{k=1}^{t}  \gamma^{k}= \frac{\gamma^{t+1}-1} {\gamma-1}-1
$$
si vede subito essendo $1-\beta^{-t} < 1$ che 
\begin{eqnarray}\label{bn}
|x|&=&\beta^e  \sum_{k=1}^{t} d_k \beta^{-k} \leq \beta^e \sum_{k=1}^{t} (\beta-1) \beta^{-k}  \nonumber \\
&=& \beta^e (\beta-1) \sum_{k=1}^{t}  \beta^{-k} = \beta^e (\beta-1)\left( \frac{1-\beta^{-(t+1)}}{1-\beta^{-1}}-1\right)  \nonumber \\
&=& \beta^e (\beta-1)\left( \frac{\beta(1-\beta^{-(t+1)})}{\beta-1}-1\right)  \nonumber \\
&=& \beta^e (\beta-1)\left( \frac{\beta(1-\beta^{-(t+1)})-(\beta-1)}{\beta-1}\right)  \nonumber \\
&=& \beta^e (1-\beta^{-t}) < \beta^e
\end{eqnarray}
e da
$$
\sum_{k=1}^{t} d_k \beta^{-k} = \frac{|x|}{\beta^e} < \frac{\beta^e}{\beta^e} < 1
$$
deduciamo che la somma relativa alla mantissa, cio\mye $\sum_{k=1}^{t} d_k \beta^{-k} $, \'e minore di $1$.
\end{itemize}.

\vspace{0.25cm}

 Se $x$ \mye il numero reale $$x=\sgn(x)  (0.d_1, \ldots d_t,\ldots)_{beta} \beta^e= \sgn(x)  \beta^e  \sum_{k=1}^{+\infty} d_k \beta^{-k}, \hspace{0.3cm} d_1 > 0$$ allora con $\myfl(x)$ denotiamo
 $$fl(x)=\sgn(x)  (0.d_1, \ldots,d_{t-1},\overline{d_t}) \beta^e= \sgn(x)  \beta^e  (\sum_{k=1}^{t-1} d_k \beta^{-k}+\overline{d_t} \beta^{-t})$$ in cui se si effettua il {\sl{troncamento}} $$\overline{d_t} = d_t$$  mentre se viene effettuato l'{\sl{arrotondamento}}.
 $$\overline{d_t} =\left \{ \begin{array}{l} d_t, \, \hspace{0.75cm} d_{t+1} < \beta/2 \\   d_t+1, \, \hspace{0.17cm}  d_{t+1} \geq  \beta/2 \end{array} \right.$$  La scelta tra troncamento o arrotondamento \mye eseguita a priori dal sistema numerico utilizzato dal calcolatore. L'ultima, \mye la pi\myu diffusa.

\vspace{0.25cm}

\begin{example}
Si consideri $$\pi=(+1) \cdot (0.314159 26535 \ldots \ldots)_{10} \cdot 10^1.$$ Se lo tronchiamo alla quarta cifra decimale otteniamo
$$fl(\pi(+1))=+1 \cdot (0.3141)_{10} \cdot 10^1,$$
mentre se lo arrotondiamo 
$$fl(\pi(+1))=+1 \cdot (0.3142)_{10} \cdot 10^1.$$
Evidentemente, per effettuare un arrotondamento, il calcolatore deve immagazzinare in un extra-bit la cifra $d_{t+1}$ (rifletterci un po').
\end{example}

\vspace{0.25cm}

Se l'esponente $e$ del numero $x$ \mye minore di $L$ si commette un'errore di {\sl{underflow}}, mentre se \mye maggiore di $U$ si commette un'errore di {\sl{overflow}}. Se invece $e \in [L,U]$, $x= \pm \beta^e  \sum_{k=1}^{+\infty} d_k \beta^{-k}$ e per qualche $t^{\ast} > t$ si ha $d_{t^{\ast}} > 0$ allora si commette un errore di {\sl{troncamento}} o {\sl{arrotondamento}} a seconda si effettui l'uno o l'altro. 

\vspace{0.25cm}

\begin{remark}
In generale la cardinalit\mya dei numeri macchina \mye molto alta, spesso ordini di grandezza maggiore di quanto sia necessaria nelle scienze applicate {\cite{TRENUM}}.

\vspace{0.25cm}

\begin{quotation}
On a computer, the interval $[1,2]$, for example, is approximated by about $10^{16}$ numbers. It is interesting to compare the fineness of this discretization with that of the discretizations of physics. In a handful of solid or liquid or a balloonful of gas, the number of atoms or molecules in a line from one point to another is on the order of $108$ (the cube root of Avogadro?s num- ber). Such a system behaves enough like a continuum to justify our definitions of physical quantities such as density, pressure, stress, strain, and temperature. Computer arithmetic, however, is more than a million times finer than this. Another comparison with physics concerns the precision to which fundamental constants are known, such as (roughly) $4$ digits for the gravitational constant $G$, $7$ digits for Planck's constant h and the elementary charge $e$, and $12$ digits for the ratio $\mu_e/\mu_{B}$ of the magnetic moment of the electron to the Bohr magneton. At present, almost nothing in physics is known to more than $12$ or $13$ digits of accuracy. Thus IEEE numbers are orders of magnitude more precise than any number in science. (Of course, purely mathematical quantities like $\pi$ are another matter.)
\end{quotation}

\end{remark}

\vspace{0.25cm}

\subsection{La precisione di macchina} Una delle costanti rilevanti all'interno del sistema floating point \mye la precisione di macchina.

\vspace{0.25cm}

\begin{definition}
Si dice {\sl{precisione di macchina}}, quel numero {\tt{eps}} che rappresenta la distanza 
tra $1$ ed il successivo numero in virgola mobile.  
\end{definition}

\vspace{0.25cm}

Calcoliamo tale valore. 
Necessariamente il numero successivo a $1$ ha per mantissa $d_1=0,\ldots,d_{t-1}=0,d_t=1$ e esponente $e=1$ (perch\mye deriva da (\ref{bn})?) e quindi \mye $$x = \beta^1 \cdot \beta^{-t}=\beta^{1-t}.$$
Si sottolinea che la precisione di macchina {\tt{eps}} non coincide in generale con il pi\myu piccolo numero positivo rappresentabile dal calcolatore, in quanto se $t \neq -L$, come usualmente accade, allora
$$
\beta^{1-t} \neq \beta^{L-1}.
$$

\vspace{0.25cm}
In Matlab ad esempio, 
\begin{itemize}
\item il valore minimo in modulo \mye $$x_{min} = 2.225073858507201 \cdot 10^{-308},$$
\item il valore massimo  \mye $$x_{max} = 1.7976931348623158 \cdot 10^{+308},$$
\item la precisione di macchina  \mye $${\mbox{\tt{eps}}}=2.220446049250313 \cdot 10^{-16}.$$
\end{itemize}
Quindi $x_{min}$ \mye molto inferiore di {\tt{eps}}.

\vspace{0.25cm}

\begin{remark}
Usualmente
\begin{itemize}
\item un numero macchina pi\myu piccolo di $x_{min}$ produce il cosidetto {\sl{underflow}},
\item un numero macchina pi\myu grande di $\pm x_{max}$ produce il cosidetto {\sl{overflow}}.
\end{itemize}
Osserviamo che esistono numeri denormalizzati, il cui range estende quello dei normalizzati. 
II più piccolo numero positivo denormalizzato è $0.494\cdot 10^{-323}$, permettendo di passare {\sl{gradualmente}} dall'underflow allo $0$.
\end{remark}

\vspace{0.25cm}

\begin{remark}
Per quanto riguarda numeri speciali
\begin{itemize}
\item per rappresentare {\tt{+inf}} si utilizza il segno pari a $0$, tutti gli esponenti uguali a $1$ e tutta la mantissa uguale a $0$;
\item per rappresentare {\tt{-inf}} si utilizza il segno pari a $1$, tutti gli esponenti uguali a $1$ e tutta la mantissa uguale a $0$;
\item per rappresentare lo $0$ si distingue $+0$ da $-0$ (generalmente utilizzati per denotare un numero troppo piccolo in modulo che viene arrotondato a $0$, pur mantenendo il segno); a parte il segno (positivo di default), tutte le altre cifre di esponente e mantissa sono uguali a $0$;
\item per rappresentare con {\tt{NaN}} il risultato di un'operazione (numerica) è stata eseguita su operandi non validi (specialmente in calcoli in virgola mobile). Secondo IEEE 754 i NaN sono rappresentati con il campo dell'esponenziale riempito di "$1$" e un numero diverso da zero nel campo della mantissa. 
\end{itemize}
\end{remark}

\vspace{0.25cm}










\subsection{Errori relativi e assoluti}

Fissato un numero da approssimare $$x^{\ast} \in {\mathbb R},$$
si definisce
{\smsp}
\begin{enumerate}
\item errore {\sl{assoluto}} tra $x$ e $x^{\ast}$ la quantit\mya
$$|x-x^{\ast}|$$
\item errore {\sl{relativo}} tra $x$ e $x^{\ast} \neq 0$ la quantit\mya 
$$ {\frac{|x-x^{\ast}|}{|x^{\ast}|}}.$$
\end{enumerate}
Si noti che se ${|x^{\ast}|}=0$, cio\mye ${x^{\ast}}=0$ per la propriet\mya del valore assoluto {\cite[p.481]{ATK}}, allora non ha senso la scrittura propria dell'errore relativo. 

\vspace{0.25cm}

\noindent
\begin{remark}
Fissato un vettore da approssimare $$x^{\ast}=(x^{\ast}_1,\ldots,x^{\ast}_n) \in {\mathbb R}^n,$$
si definisce
{\smsp}
\begin{enumerate}
\item errore {\sl{assoluto}} tra $x$ e $x^{\ast}$ la quantit\mya
$$\|x-x^{\ast}\|$$
dove $$\|y\|=\sqrt{\sum_{i=1}^n y_i^2}, \, \, y=(y_1,\ldots,y_n);$$
{\smsp}
\item errore {\sl{relativo}} tra $x$ e $x^{\ast} \neq 0$ la quantit\mya 
$$ {\frac{\|x-x^{\ast}\|}{\|x^{\ast}\|}}.$$
\end{enumerate}
Si noti che se ${\|x^{\ast}\|}=0$, cio\mye ${x^{\ast}}=0$ per la propriet\mya delle norme {\cite[p.481]{ATK}}, allora non ha senso la scrittura propria dell'errore relativo. Inoltre se il vettore ha un solo componente, allora la norma $\|\cdot\|$ coincide con l'usuale valore assoluto $|\cdot|$. $\qed$
\end{remark}

\vspace{0.25cm}

Sia $$x = \sgn(x) \beta^e  \sum_{k=1}^{+\infty} d_k \beta^{-k} \neq 0, \hspace{0.3cm} d_1 >0$$ un numero reale e $$fl(x)=\sgn(x)  (0.d_1, \ldots,d_{t-1},\overline{d_t}) \beta^e= \sgn(x)  \beta^e  (\sum_{k=1}^{t-1} d_k \beta^{-k}+\overline{d_t} \beta^{-t}), \hspace{0.3cm} d_1 >0$$ la sua approssimazione in virgola mobile.

L'errore relativo compiuto nell'approssimare $x$ con $fl(x)$, visto che $|x| \geq |{d_1} \beta^{e-1}| \geq |\beta^{e-1}|$ poich\mye $d_1 >0$, verifica
\begin{eqnarray}
\frac{|x-fl(x)|}{|x|} &\leq&  \frac{| \sgn(x) \beta^e  \sum_{k=1}^{+\infty} d_k \beta^{-k}-(\sgn(x)\beta^e  (\sum_{k=1}^{t-1} d_k \beta^{-k}+\overline{d_t} \beta^{-t}))|}{ |\beta^{e-1}| } \nonumber \\
&=& \beta \cdot \left| \sum_{k=t}^{+\infty} d_k \beta^{-k}-\overline{d_t} \beta^{-t} \right| 
\end{eqnarray}


Nel caso si effettui il troncamento abbiamo ${d_t} = \overline{d_t}$ e quindi, essendo ${\mbox{\tt{eps}}}=\beta^{1-t}$ la precisione di macchina
$$  \beta \cdot \left| \sum_{k=t}^{+\infty} d_k \beta^{-k}-\overline{d_t} \beta^{-t} \right|  =  \beta \cdot \left| \sum_{k=t+1}^{+\infty} d_k \beta^{-k}\right| < \beta \cdot(\beta \cdot \beta^{-t-1})=\beta^{1-t}={\mbox{\tt{eps}}}$$
mentre se si effettua l'arrotondamento abbiamo dopo alcuni calcoli che 
$$\frac{|x-fl(x)|}{|x|} < \frac{\beta^{1-t}}{2}= \frac{{\mbox{\tt{eps}}}}{2}.$$

Il minor errore relativo compiuto nell'approssimare $x$ con $fl(x)$ via arrotondamento, suggerisce perch\mye quest'ultimo sistema sia pi\myu utilizzato.

\vspace{0.25cm}

\begin{remark}
Da questo si deduce che i 
numeri in virgola mobile non sono dunque un insieme denso, come sono 
invece i numeri reali, bens\myi un insieme discreto di valori sull'asse 
reale, in posizioni non equispaziate (se cos\myi fosse, visto che i numeri pi\myu grandi sono interi, pure {\tt{eps}} lo sarebbe). 
\end{remark}








\section{Operazioni con i numeri macchina}
Per quanto concerne le quattro operazioni aritmetiche, si pu\myo dimostrare che la divisione e la moltiplicazione di tipo floating-point producono un errore relativo piccolo rispetto al risultato esatto, cosa che non succede sempre per somma e sottrazione.

\vspace{0.25cm}

\noindent
Pi\myu precisamente se ${\mbox{fl}}(x)$ \mye il numero macchina che {\sl{corrisponde}} a $x$, denotiamo con 

\begin{eqnarray}
x \oplus y &:=& \myfl (\myfl(x)+\myfl(y)) \\
x \ominus y &:=& \myfl (\myfl(x)+\myfl(y)) \\
x \otimes y &:=& \myfl (\myfl(x) \cdot \myfl(y)) \\
x \oslash y &:=& \myfl (\myfl(x):\myfl(y)) 
\end{eqnarray}

\noindent
e per ${\odot}$ una delle operazioni precedentemente introdotte, cio\mye una tra $\oplus$, $\ominus$, $\otimes$, $\oslash$ corrispondenti a {\tt{op}} cio\mye una tra $+$, $-$, $\cdot$, $:$. Inoltre porremo
\begin{eqnarray}
\epsilon_x:=\frac{|x-\myfl(x)|}{|x|} \\
\epsilon_{x,y}^{\odot}:=\frac{| (x \, {\mbox{op}}  \, y )-(x \, {\odot}  \, y )|}{|x \, {\mbox{op}} \, y|}
\end{eqnarray}

\begin{remark}
Osserviamo, che a volte nei manuali si trova qualcosa di diverso. Ad esempio, 
$$
x \oplus y = \myfl (x+y) 
$$
In realt\mya il sistema floating point, utilizzando ulteriori bits nei calcoli, garantisce che sia 
$$
x \oplus y = \myfl (\myfl(x)+\myfl(y)) = \myfl (x+y).
$$ 
In effetti, da {\cite{GOLD}}

\vspace{0.25cm}

\begin{quotation}
The IEEE standard requires that the result of addition, subtraction, multiplication and division be exactly rounded. That is, the result must be computed exactly and then rounded to the nearest floating-point number (using round to even). The section Guard Digits pointed out that computing the exact difference or sum of two floating-point numbers can be very expensive when their exponents are substantially different. That section introduced guard digits, which provide a practical way of computing differences while guaranteeing that the relative error is small. However, computing with a single guard digit will not always give the same answer as computing the exact result and then rounding. By introducing a second guard digit and a third sticky bit, differences can be computed at only a little more cost than with a single guard digit, but the result is the same as if the difference were computed exactly and then rounded . Thus the standard can be implemented efficiently.
\end{quotation}

\vspace{0.25cm}

\noindent Altre volte si trova equivalentemente {\cite{TRENUM}}

\vspace{0.25cm}

\begin{quotation}
if $*$ is one of these four operations in its ideal form and $\odot$ is the same operation as realized on the computer, then for any floating-point numbers $x$ and $y$ , assuming that there is no underflow or overflow,
$$x \odot y=(x*y)(1+\epsilon).$$
Here $\epsilon$ is a very small quantity, no greater in absolute value than a number known as machine epsilon, denoted by $\epsilon_{mach}$, that measures the accuracy of the computer. In the IEEE system, 
$\epsilon_{mach} = 2^{-53} \approx 1.1 \times 10^{-16}$.
\end{quotation}
\end{remark}

\vspace{0.25cm}

\subsection{Propriet\'a commutativa, associativa e distributiva delle operazioni floating point} Per prima cosa osserviamo che alcune propriet\mya caratteristiche di $+$, $-$, $\cdot$ e $:$ non valgono per i corrispettivi $\oplus$, $\ominus$, $\otimes$, $\oslash$. Infatti per la somma $\oplus$ e il prodotto $\otimes$ vale la propriet\mya commutativa ma non vale quella associativa, e nemmeno quella distributiva.

\vspace{0.25cm}

\begin{example}
Mostriamo come non valga la propriet\mya associativa della somma. Consideriamo $F(10,4,-10,+10)$ con {\sl{fl}} basato sul troncamento. Sia $a=2000$, $b=2.5$, $c=7.8$. Allora:
\begin{eqnarray}
(a \oplus b) \oplus c &=& 0.2002 \cdot 10^4 \oplus 0.7800 \cdot 10^1 = 0.2009 \cdot 10^4 \nonumber \\
 a \oplus (b \oplus c)&=& 0.2000 \cdot 10^4 \oplus 0.1030 \cdot 10^2 = 0.2010 \cdot 10^4 \nonumber
\end{eqnarray}
e quindi $(a \oplus b) \oplus c \neq  a \oplus (b \oplus c)$.

{\vspace{0.25cm}}

Per quanto riguarda la propriet\mya distributiva, in $F(10,4,-10,+10)$, posti $a=0.800000 \cdot 10^9$, $b=0.500009 \cdot 10^4$, $c=0.500008 \cdot 10^4$
\begin{eqnarray}
a \otimes (b \ominus c) &=& 0.800000 \cdot 10^9 \otimes 0.1 \cdot 10^{-2} = 0.800000 \cdot 10^9 \nonumber \\
(a \otimes b) \ominus (a \otimes c) &=& (0.800000 \cdot 10^9 \otimes b=0.500009 \cdot 10^4) \ominus \nonumber \\ &&\ominus (0.800000 \cdot 10^9 \otimes b=0.500008 \cdot 10^4) \rightarrow {\mbox{overflow}}
\end{eqnarray}
(l'overflow \mye dovuto alla sequenza di operazioni da compiere).
\end{example}

\subsection{Errori nelle operazioni e loro propagazione}
Consideriamo ora come si propagano gli errori in $\oplus$, $\ominus$, $\otimes$, $\oslash$. Si mostra che (cf. {\cite[p.78]{COM}})

\vspace{-0.5cm}

\begin{eqnarray}
\epsilon_{x,y}^{\oplus} &\leq&  \left| {\frac{x}{x+y} } \right| \epsilon_x +  \left| {\frac{y}{x+y} } \right| \epsilon_y \\
\epsilon_{x,y}^{\ominus} &\leq&   \left| {\frac{x}{x-y} } \right| \epsilon_x +  \left| {\frac{y}{x-y} } \right| \epsilon_y \\
\epsilon_{x,y}^{\otimes} &\lessapprox& \epsilon_x + \epsilon_y \\
\epsilon_{x,y}^{\oslash} &\leq&  \left| \epsilon_x - \epsilon_y \right|
\end{eqnarray}


\vspace{0.25cm}

Proviamo di seguito sotto opportune ipotesi il caso della somma, sottrazione e della moltiplicazione, lasciando al lettore quello della divisione. Con $a \lessapprox b$ si intende $a$ {\sl{inferiore}} o uguale circa a $b$.

{\vspace{0.25cm}}

\subsubsection{Errore nella somma}
Per quanto riguardo la somma, supposti $x+y \neq 0$, $x \neq 0$, $y \neq 0$, e $\myfl (\myfl(x)+\myfl(y))= \myfl(x)+\myfl(y)$ (cio\mye $\myfl(x)+\myfl(y)$ numero macchina, richiesta non eccessiva, vuol dire semplicemente che il risultato \mye ancora nell'insieme dei numeri macchina), ricordato che $\epsilon_x=\frac{|x-\myfl(x)|}{|x|}$ed $\epsilon_y=\frac{|y-\myfl(y)|}{|y|}$, abbiamo
\begin{eqnarray}
\epsilon_{x,y}^{\oplus}&=&\frac{|(x+y)-(x \oplus y)|}{|x+y|}=\frac{|(x+y)-\myfl (\myfl(x)+\myfl(y))|}{|x+y|}=\frac{|(x-\myfl(x))+(y-\myfl(y))|}{|x+y|} \nonumber \\
&\leq& \frac{|x|}{|x+y|}  \cdot  \frac{|x-\myfl(x)|}{|x|}+  \frac{|y|}{|x+y|} \cdot \frac{|y-\myfl(y)|}{|y|}   =   \left| {\frac{x}{x+y} } \right| \epsilon_x +  \left| {\frac{y}{x+y} } \right| \epsilon_y
\end{eqnarray}
Il caso in cui $\myfl (\myfl(x)+\myfl(y)) \neq \myfl(x)+\myfl(y)$ \mye pi\myu complicato ma i risultati restano simili dal punto di vista teorico.

Il punto chiave \mye che per $x+y \approx 0$ abbiamo che per piccoli errori sui dati $\epsilon_x$, $\epsilon_y$ possono aversi errori rilevanti in $\epsilon_{x,y}^{\oplus}$. 

\vspace{0.25cm}

Per quanto riguarda la sottrazione, con una analisi simile, si vede che
\begin{eqnarray}
\epsilon_{x,y}^{\ominus}&=&\frac{|(x-y)-(x \ominus y)|}{|x-y|}=\frac{|(x-y)-\myfl (\myfl(x)-\myfl(y))|}{|x-y|}=\frac{|(x-\myfl(x))-(y-\myfl(y))|}{|x-y|} \nonumber \\
&\leq& \frac{|x|}{|x-y|}  \cdot  \frac{|x-\myfl(x)|}{|x|}+  \frac{|y|}{|x-y|} \cdot \frac{|y-\myfl(y)|}{|y|}   =  \left| {\frac{x}{x-y} } \right| \epsilon_x +  \left| {\frac{y}{x-y} } \right| \epsilon_y
\end{eqnarray}

Nuovamente, per $x-y \approx 0$ abbiamo che per piccoli errori sui dati $\epsilon_x$, $\epsilon_y$ possono aversi errori rilevanti in $\epsilon_{x,y}^{\ominus}$. Qualora ci\myo accada, si parla di {\sl{fenomeno di cancellazione}}.

\vspace{0.25cm}

\begin{example}
Consideriamo il sistema $F(10,5,-5,+5)$ (con arrotondamento) e i due numeri 
$$a=0.73415776, {\hspace{0.5cm}}b=0.73402350.$$
Naturalmente, 
$$\myfl{(a)}=0.73416, {\hspace{0.5cm}}\myfl{(b)}=0.73402$$
e quindi 
$$a-b=0.00013426, {\hspace{0.5cm}} a \ominus b=0.00014=10^{-3} \cdot 0.14000$$
con un errore relativo pari a 
$$|(a-b)- (a \ominus b)|/|a-b|\approx 4.3/100$$
cio\mye molto grande (oltre il 4 per cento!).
\end{example}



\subsubsection{Errore nel prodotto}

Nel caso del prodotto,  
\begin{itemize}
\item supposti $x \cdot y \neq 0$, $x \neq 0$, $y \neq 0$, $|fl(y)|/y \approx 1$,
\item supposto $\myfl (\myfl(x) \cdot \myfl(y))= \myfl(x) \cdot \myfl(y)$ (cio\mye $\myfl(x) \cdot \myfl(y)$ numero macchina), 
\item ricordato che $\epsilon_x=\frac{|x-\myfl(x)|}{|x|}$ed $\epsilon_y=\frac{|y-\myfl(y)|}{|y|}$, 
\end{itemize}
abbiamo
\begin{eqnarray}
\epsilon_{x,y}^{\otimes}&=&\frac{|x\cdot y-x \otimes y|}{|x \cdot y|}=\frac{|x\cdot y-\myfl (\myfl(x) \cdot \myfl(y))|}{|x\cdot y|}=\frac{|x\cdot y - \myfl(x) \cdot \myfl(y)|}{|x\cdot y|} \nonumber \\
&=& \frac{|x\cdot y - x\cdot \myfl (y)+ x\cdot \myfl (y)- \myfl(x) \cdot \myfl(y)|}{|x\cdot y|} \nonumber \\
&\leq&  \frac{| x\cdot (y-\myfl (y))|+|\myfl(y) \cdot (x-\myfl(x)) | }{|x\cdot y|}  \nonumber \\
&=& \frac{| y-\myfl (y)|}{| y|} + \frac{|\myfl(y) \cdot (x-\myfl(x)) | }{|x\cdot y|} \nonumber \\
& \approx &  \epsilon_x + \epsilon_y
\end{eqnarray}

\noindent 
dove si ricorda che $$ \frac{|\myfl(y)| }{|x\cdot y|}  \approx  \frac{1}{|x|} .$$

\vspace{0.25cm}

\begin{remark}
Per stabilit\mya di una generica operazione macchina ${\mbox{op}}$ (si pensi ad esempio $*$ come alla somma oppure al prodotto), si intende che se $\odot$ \mye l'omologa versione floating-point allora esiste $C > 0$ tale che
$$
\epsilon_{x,y}^{\odot} \leq C (\epsilon_{y} {\mbox{ op }} \epsilon_{y})
$$
per ogni $x,y$ per cui la scrittura abbia significato.

Ne consegue che la somma e la sottrazione non sono stabili mentre lo sono il prodotto e la divisione \cite[p.10]{QUAold}.
\end{remark}

\vspace{0.25cm}

\begin{remark}
Osserviamo che dalla formula di Taylor in due variabili, posto per brevit\mya di notazione $f^{\prime}_x=\partial f / \partial x$, $f^{\prime}_y=\partial f / \partial y$,
\begin{equation}\label{T2D}
f(x,y) \approx f({\overline{x}},{\overline{y}})+f_x^{\prime}({\overline{x}},{\overline{y}}) (x-\overline{x})+f_y^{\prime}({\overline{x}},{\overline{y}}) (y-\overline{y}) 
\end{equation}
necessariamente, se $x,y,f(x,y) \neq 0$ allora posti $\epsilon_x=|x-\overline{x}|/|x|$, $\epsilon_y=|y-\overline{y}|/|y|$ gli errori relativi sui dati
\begin{eqnarray}\label{T2D1}
\frac
{|f(x,y)- f( {\overline{x}} , {\overline{y}})|}
{|f(x,y)|}  
&\lessapprox& 
\frac{|f_x^{\prime}({\overline{x}},{\overline{y}})||x-\overline{x}||x|}
{|f(x,y)||x|}+\frac{ |f_y^{\prime}({\overline{x}},{\overline{y}})||y-\overline{y}||y|}{|f(x,y)||y|} \nonumber \\
&=&\frac{|f_x^{\prime}({\overline{x}},{\overline{y}})||x|}
{|f(x,y)|} \cdot \epsilon_x+\frac{ |f_y^{\prime}({\overline{x}},{\overline{y}})||y|}{|f(x,y)|} \cdot \epsilon_y
\end{eqnarray}
Da (\ref{T2D1}), si possono ottenere risultati simili ai precedenti per $\oplus$, $\otimes$. Uno dei vantaggi \mye che si possono studiare operazioni pi\myu generali come ad esempio cosa succeda se si esegue $exp(x+y)$.
\end{remark}

\vspace{0.25cm}




\section{Alcune problematiche numeriche}
Si consideri un problema scientifico che abbia a che fare con una procedura di calcolo. Sussitono varie tipologie di errore nei calcoli:
\begin{itemize}
\item errori di modelizzazione: vengono utilizzate equazioni per rappresentare un evento fisico, e il modello matematico \mye solo una semplificazione di quello reale;
\item errori di tabulazione: vengono immessi dal programmatore o dal programma dati errati (cosa che ad esempio era comune una volta quando bisognava tabulare manualmente numeri con molte cifre); un errore di questo tipo ha causato nel 1996 la perdita del razzo Ariane 5 per una cattiva conversione di un numero;
\item errori dovuti a misure fisiche: sono dovuti a una imprecisa osservazione sperimentale, spesso dovuta agli strumenti utilizzati, ad esempio la velocit\'a della luce nel vuoto \mye
$$
c=(2.997925+\epsilon)\cdot 10^{10} cm/s, \,\,\, |\epsilon| \leq 0.000003
$$
e quindi nei calcoli che hanno a che vedere con $c$ necessariamente si compie un errore;
\item errori di rappresentazione e di aritmetica di macchina: sono errori dovuti all'arrotondamento o troncamentodi un numero, e sono inevitabili quando si usa l'aritmetica floating-point; sono la sorgente principale di errore di alcuni problemi come ad esempio la soluzione di sistemi lineari. 
\end{itemize}

Di seguito vediamo alcune problematiche dovute agli errori di rappresentazione e di aritmetica di macchina. Per una buona serie di esempi si veda {\cite[p.46]{ATKHAN}}.


\section{Valutazione di una funzione}
Nel valutare una funzione continua $f: \Omega \subseteq \mathbb{R}  \rightarrow \mathbb{R} $ in un punto $x \in \Omega$ bisogna tener conto che il dato $x$ pu\myo essere per varie ragioni affetto da errori, tipicamente di misura o di rappresentazione al calcolatore, cos\myi da essere approssimato con un certo ${\overline{x}}$. Di conseguenza, invece di $f(x)$ si valuta $f({\overline{x}})$ e ci si domanda se ci siano considerevoli differenze.

Una funzione $f$ risulta difficile da valutare al calcolatore nel punto $x \neq 0$ in cui $f(x) \neq 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)|.$$ Quindi \mye importante discutere la quantit\mya
$$ {\cal{K}}(f,x,x_c)=\frac{ |f(x)-f(x_c)|/|f(x)|}{|x-x_c|/|x|}.$$
La valutazione di $f$ in un punto $x$ si dice {\sl{bencondizionata}} se a piccole variazioni relative sui dati corrispondono piccole variazioni sui risultati,  {\sl{malcondizionata}}  se ci\myo non avviene.
Se $f$ \mye derivabile con continuit\mya nel pi\myu piccolo intervallo $\cal{I}$ aperto contenente $x$ ed $x_c$, per il teorema della media, essendo $$f(x)-f(x_c)=f^{\prime}(\xi) \cdot (x-x_c), \,\, \xi \in {\cal I},$$ e quindi da $|f(x)-f(x_c)|/|x-x_c| = f^{\prime}(\xi) \approx f^{\prime}(x)$ ricaviamo facilmente che
$$ {\cal{K}}(f,x,x_c) \approx {\cal{K}}(f,x) := \frac{|  x \cdot f^{\prime}(x)|}{|f(x)|} .$$
La quantit\mya ${\cal{K}}(f,x)$ si chiama {\sl{condizionamento}} di $f$ nel punto $x$. Pi\myu piccola risulta ${\cal{K}}(f,x)$ e meno la funzione amplifica un piccolo errore sul dato $x$. Naturalmente la funzione pu\myo essere {\sl{bencondizionata}} su un certo dato $x_1 \in \Omega$ e {\sl{malcondizionata}} su un certo altro $x_2  \in \Omega$.
Osserviamo che il condizionamento di una funzione non dipende dall'algoritmo con cui viene valutata, ma \mye inerente alla funzione stessa $f$ e al dato $x$.

\vspace{0.25cm}

\begin{example}
Data la funzione
$$f_1(x)=1-\sqrt{1-x^2}$$
calcoliamo analiticamente il {\sl{condizionamento}}
$${\cal{K}}(f_1,x)=\frac{|  x \cdot f_1^{\prime}(x)|}{|f_1(x)|} $$
Ricordiamo che
$$
f_1(x)={\frac{1}{2 \sqrt{1-x^2}}} (-2x) 
$$

\noindent da cui per $x \in (-1,1)$
$${\cal{K}}(f_1,x)=\frac{|  x \cdot f_1^{\prime}(x)|}{|f_1(x)|} = \frac{|  x \cdot{\frac{1}{2 \sqrt{1-x^2}}} (-2x) |}{|1-\sqrt{1-x^2}|} = \frac{x^2}{|1-\sqrt{1-x^2}|\sqrt{1-x^2} }.$$

\noindent Dal grafico di tale funzione nella Figura \ref{fig1_cond}, si vede che la funzione $f_1$ \mye ben condizionata in $(-0.9,0.9)$, ma non lo \mye per valori prossimi a $1$.

 \begin{figure}[h!]\label{fig1_cond}
 \centering
   {\includegraphics[scale=0.15,clip]{cond_1C.jpg}}
   {\includegraphics[scale=0.15,clip]{cond_1_zoomC.jpg}}
 \caption{A sinistra il grafico di ${\cal{K}}(f_1,x)$ in $(-1,1)$, a destra quello di ${\cal{K}}(f_1,x)$ in $(1-10^{-4},1)$ (malcondizionamento per $|x| \approx 1$).}
 \end{figure}
 \end{example}

\vspace{0.25cm}

 \begin{example}
 Consideriamo quale esempio
 $$f_2(x)=1-x.$$
 Non \mye difficile vedere che $f^{\prime}_2(x)=-1$ e quindi
 $${\cal{K}}(f_2,x)=\frac{|  x \cdot f_2^{\prime}(x)|}{|1-x|}=\frac{|  x \cdot (-1)|}{|1-x|}=\frac{|  x \cdot (-1)|}{|1-x|}.$$
 \noindent Dal grafico di tale funzione nella Figura \ref{fig2_cond}, si vede che la funzione $f_2$ \mye ben condizionata ovunque eccetto per valori prossimi a $1$. Si osservi che per vicini a tale valore, in effetti, si ha il fenomeno di cancellazione.

 \begin{figure}[h!]\label{fig2_cond}
 \centering
   {\includegraphics[scale=0.15,clip]{cond_2.jpg}}
 \caption{Il grafico di ${\cal{K}}(f_2,x)$ in $(-2,2)$ (malcondizionamento per $x \approx 1$).}
 \end{figure}
 
\end{example}

\vspace{0.25cm}


\section{Stabilit\myas}
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 {\sl{stabile}} se amplifica {\sl{poco}} gli errori di arrotondamento introdotti dalle singole operazioni. 
Il problema non \mye necessariamente la valutazione di una funzione continua in un punto, ma pu\myo essere la risoluzione di una equazione di secondo grado, il calcolo di $\pi$ 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 \mye affrontato con un algoritmo instabile. 

\vspace{0.25cm}

\subsection{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\mya della soluzione e come rimediare.
Dato il polinomio di secondo grado $x^2 + 2\,p x -q$, con 
${\sqrt{p^2+q}}\geq 0$ calcoliamo la 
radice 
\begin{equation}\label{quadeq}
y=-p+{\sqrt{p^2+q}}.
\end{equation}

\noindent
Osserviamo che essendo ${\sqrt{p^2+q}}\geq 0$ le radici 
\begin{equation}\label{quadeq2}
y=-p \pm {\sqrt{p^2+q}}.
\end{equation}
dell'equazione sono reali. La soluzione descritta in (\ref{quadeq}) \mye la maggiore delle 2, ed \mye non negativa se e solo se $q \geq 0$. Consideriamo la funzione che fissato $q \geq 0$, mappa un numero reale $p$ in $-p \pm {\sqrt{p^2+q}}$.
Si osserva subito che (\ref{quadeq}) \mye potenzialmente instabile per $p \gg 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 ({\ref{quadeq}}). 
In altri termini
\begin{equation}\label{quadeqr}
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}})} }
\end{equation}

\noindent
Ricordiamo ora che un problema si dice {\sl{bencondizionato}} (o {\sl{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 {\sl{algoritmo bencondizionato}} (o {\sl{algoritmo malcondizionato}}) anche se \mye pi\myu usuale il termine di stabilit\mya 
{\cite[p. 66]{COM}}.

Seguendo 
{\cite[p. 10]{COM}}, {\cite[p. 78]{COM}}, il problema (e non l'algoritmo!) \mye {\sl{bencondizionato}}  per $q>0$ e {\sl{malcondizionato}} per $q \approx -p^2$. 

\noindent
Usando dei classici ragionamenti dell'analisi numerica si mostra che (cf. {\cite{STO}}, p. 21, {\cite{COMESE}}, p. 11)
{\smsp}
\begin{enumerate}
\item il primo algoritmo (\ref{quadeq2}) non \mye {\sl{numericamente stabile}} qualora $p \gg q > 0$; 
{\smsp}
\item  il secondo algoritmo (\ref{quadeqr}) \mye {\sl{numericamente stabile}} qualora $p \gg q  > 0$. 
\end{enumerate}

\vspace{0.25cm}

\noindent Vediamo alcuni casi rilevanti, effettuati in doppia precisione.

\vspace{0.25cm}

\begin{itemize}
\item In {\cite[p.22]{STO}}, si suggerisce un test interessante per $$p=1000, \, q=0.018000000081$$ la cui soluzione esatta \mye $0.9 \,\cdot\,10^{-5}$. Si ha $p \gg q$. Lavorando in doppia precisione otteniamo

\vspace{0.25cm}

\begin{tabular}{| c | c | c |}
\hline
{\mbox{algoritmo}}  & {\mbox{valore}} & {\mbox{err. rel.}} \\
\hline
$1$ &  $0.0000089999999772772$ & $2.52\cdot 10^{-9}$ \\
$2$ & $ 0.0000090000000000000 $ & $0$ \\
\hline
\end{tabular}

\vspace{0.25cm}

\item Secondo {\cite[p.11]{COMESE}}, \mye notevole l'esperimento in cui $$p=4.999999999995 \cdot 10^{+4}, \, q=10^{-2}$$ avente soluzione esatta $10^{-7}$.
Nuovamente $p \gg q$.

\vspace{0.25cm}

\begin{tabular}{| c | c | c |}
\hline
{\mbox{algoritmo}}  & {\mbox{valore}} & {\mbox{err. rel.}} \\
\hline
$1$ &  $0.0000001000007614493$ & $7.61\cdot 10^{-6}$ \\
$2$ & $ 0.0000001000000000000 $ & $0$ \\
\hline
\end{tabular}

\vspace{0.25cm}

\end{itemize}



\subsection{Esempio 2: Approssimazione di $\pi$}

Storicamente sono state scoperte diverse successioni convergenti sempre pi\myu rapidamente a $\pi$ (cf. {\cite{WIKPI}}).
In questa sezione ci interesseremo a tre di queste mostrando sia come le velocit\mya di convergenza possano essere diverse, sia come possano insorgere questioni di instabilit\mya {\cite[p. 16]{FRO}}.

\vspace{0.25cm}

\noindent
Studiamo le successioni $\{u_n\}$, $\{z_n\}$, definite rispettivamente come
$$\left\{ \begin{array} {l}  
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{array} \right.$$
e
$$\left\{ \begin{array} {l}  
z_1=1, \, \, z_2=2 \\
z_{n+1} = 2^{n- { \frac {1} {2} } }      { \sqrt { 1 - \sqrt{1-4^{1-n} 
\cdot z_n^2}   } } 
\end{array} \right.$$
che {\sl{teoricamente}} convergono a $\pi=3,14159 26535 89793 23846 26433 83279 50288 41971 69399 \ldots$.

Implementiamo di seguito una terza successione, diciamo $\{y_{n}\}$, che si ottiene 
{\sl{razionalizzando}} $z_n$, cio\mye moltiplicando numeratore e denominatore 
per 
$$ { \sqrt { 1 + \sqrt{1-4^{1-n} \cdot z_n^2}   } } $$
e calcoliamo $u_m$, $z_m$ e $y_m$ per $m=2,3,\ldots,40$ (che teoricamente dovrebbero approssimare $\pi$).

Si ottengono i seguenti risultati, che tabuliamo per semplicit\mya esclusivamente per $n=5,10,15,20,25,30,35,40$.

\vspace{0.35cm}

\begin{tabular}{| c | c | c | c |}
\hline
$m$  & $u_n$ & $z_n$ & $y_n$ \\
\hline
 	 $  5$ & $2.963387701038571$ & $3.121445152258050$ & $3.121445152258053$ \\
 	 $ 10$ & $3.049361635982070$ & $3.141572940365566$ & $3.141572940367093$ \\
 	 $ 15$ & $3.079389826032086$ & $3.141592633463248$ & $3.141592634338565$ \\
 	 $ 20$ & $3.094669524113704$ & $3.141594125195191$ & $3.141592653570997$ \\
 	 $ 25$ & $3.103923391700576$ & $3.142451272494133$ & $3.141592653589779$ \\
 	 $ 30$ & $3.110128728141262$ & $4.000000000000000$ & $3.141592653589798$ \\
 	 $ 35$ & $3.114578862293132$ & $0.000000000000000$ & $3.141592653589798$\\
 	 $ 40$ & $3.117926198299378$ & $0.000000000000000$ & $3.141592653589798$\\
\hline
\end{tabular}

\vspace{0.35cm}

 \begin{figure}[h!]\label{pigreco}
 \centering
   {\includegraphics[scale=0.35,clip]{pigreco.eps}}
 \caption{Le tre successioni  $u_n$ (puntino nero), $z_n$ (punto magenta), $y_n$ (cerchietto rosso). Risulta evidente che $z_n$ diverge, nonostante fino a $n \approx 14$ dia circa gli stessi valori di $y_n$ che invece converge.}
 \end{figure}
 
\vspace{0.35cm}

Dalla figura e dalla tabella, consegue che
\begin{enumerate}
\item la convergenza della prima successione \mye estremamente lenta;
\item la seconda successione non convergente numericamente, seppure converga teoricamente; il principale problema \mye la cancellazione che si ha nell'argomento della radice quadrata pi\myu esterna di 
$$
z_{n+1} = 2^{n- { \frac {1} {2} } }      { \sqrt { 1 - \sqrt{1-4^{1-n} \cdot z_n^2}   } } ;
$$
A tal proposito, nella tabella che segue, evidenziamo i valori assunti da $\sqrt{1-4^{1-n} \cdot z_n^2}$, per $n=5,10,15,\ldots,40$.

\vspace{0.35cm}

\begin{center}
\begin{tabular}{| c | c |} %\label{tabrad}
\hline
$n$  & $     \sqrt{1-4^{1-n} \cdot z_n^2} $ \\
\hline
 	 $ 5$ & $9.8078528040323054160154470e-01$ \\
 	 $10$ & $9.9998117528260110908888691e-01$ \\
 	 $15$ & $9.9999998161642922323011362e-01$ \\
 	 $20$ & $9.9999999998204724960260137e-01$ \\
 	 $25$ & $9.9999999999998245847621092e-01$ \\
 	 $30$ & $1.0000000000000000000000000e+00$ \\
 	 $35$ & $1.0000000000000000000000000e+00$ \\
 	 $40$ & $1.0000000000000000000000000e+00$ \\
\hline
\end{tabular} 	 
\end{center}
\vspace{0.35cm}

\item la terza successione risulta rapidamente convergente.
\end{enumerate}

\subsection{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 {\cite[Esercizio 1.9, p.35]{QUA}}. Problemi simili con una precisa discussione della propagazione dell'errore sono trattati pure in {\cite[p. 23, problema 11]{FRO}}

Sia 
\begin{equation}\label{321}
I_n=e^{-1} \int_0^1 x^n \, e^x \, dx
\end{equation}

\noindent
Per $n=0$ si ha $$I_0=e^{-1} \int_0^1 \, e^x \, dx=e^{-1}(e^{1}-1).$$
Per $n=1$, \mye facile verificare che, essendo
$$\int x \, \exp(x) dx=(x-1)\, \exp(x)+c,$$
dal secondo teorema del calcolo integrale (cf. {\cite{WIKCI}})
si ha che $I_1=e^{-1}$ e pi\myu in generale integrando per parti che
\begin{equation}\label{322}
I_{n+1}=e^{-1} \left( x^{n+1} \, e^x \mid_0^1 - (n+1) \int_0^1 x^n \, e^x 
\, dx  \right) =1-(n+1)\, I_n. 
\end{equation}

\noindent
Da (\ref{321}) essendo l'integranda $x^n \exp{x} >0 $ per $x \in (0,1]$
si ha $$I_n >0.$$ Inoltre $x \leq 1$ implica $x^2 \leq x$ e pi\myu in generale $x^{n+1} \leq x^n$ da cui $x^{n+1} \exp(x) \leq x^n \exp(x)$ e quindi $$I_{n+1} \leq I_n.$$
La successione $I_n$ \mye positiva e non crescente e quindi ammette limite $L$ finito.
Da ({\ref{322}}), calcolando il limite $L$ per $n \rightarrow \infty$ 
ad ambo i membri si ottiene portando $1$ a primo membro
$$L-1=-\lim_n  (n+1)\, I_n.$$
L'unica possibilit\mya affinch\mye $\lim_n  (n+1)\, I_n$ esista e sia finito \mye 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=+\infty,$$ il che \mye assurdo essendo $L < + \infty$ e quindi $L-1 < +\infty$.
{\smsp}

\begin{enumerate}
\item Come primo tentativo calcoliamo $I_n$ per $n=1,\ldots, 100$ mediante la successione {\sl{in avanti}}
$$\left\{ \begin{array}{l}  
s_1=e^{-1}\\
s_{n+1} = 1-(n+1) \, s_{n}
\end{array} \right.$$

\item Come alternativa, fissato $m=500$, si calcoli la successione {\sl{all'indietro}}
$\{t_n\}_{n=1,\ldots,100}$ definita come 
$$\left\{ \begin{array}{l}  
t_{2m}=0\\
t_{n-1} = {\frac {1-t_n} {n}}
\end{array} \right.$$
Si osservi che per raggiungere tale obiettivo bisogna calcolare i termini 
$$t_m, \,t_{m-1},\ldots,t_{100},t_{99},\ldots,t_2,t_1.$$

Notiamo che tale successione all'indietro deriva dall'osservare che per $n$ sufficientemente grande $I_n \approx 0$. Quindi, posto per $n$ sufficientemente grande $I_n=0$ riscrivendo la successione in avanti come successione all'indietro (cio\mye $I_n$ in funzione di $I_{n+1}$), abbiamo 
$$I_{n} = {\frac {1-I_{n+1}} {n+1}}.$$
\end{enumerate}

Dopo aver eseguito i calcoli otteniamo per $n=10,20,\ldots,100$ la seguente tabella

\vspace{0.35cm}

\begin{center}
\begin{tabular}{| c | c | c |} %\label{tab45}
\hline
$n$  & $     \sqrt{1-4^{1-n} \cdot z_n^2} $ \\
\hline
 	 $ 10$ & $8.387707005829270e-02$ & $8.387707010339417e-02$ \\
 	 $ 20$ & $-3.019239488558378e+01$ & $4.554488407581805e-02$\\
 	 $ 30$ & $-3.296762455608386e+15$ & $3.127967393216808e-02$\\
 	 $ 40$ & $-1.014081007335147e+31$ & $2.382272866903348e-02$\\
 	 $ 50$ & $-3.780092725853228e+47$ & $1.923775443433938e-02$\\
 	 $ 60$ & $-1.034194991387092e+65$ & $1.613316466384628e-02$\\
 	 $ 70$ & $-1.488787166378969e+83$ & $1.389153285400528e-02$\\
 	 $ 80$ & $-8.895191521232202e+101$ & $1.219691454831806e-02$\\
 	 $ 90$ & $-1.846559775027300e+121$ & $1.087083605943215e-02$\\
 	 $100$ & $-1.159928542966359e+141$ & $9.804855005376052e-03$\\
\hline
\end{tabular} 	 
\end{center}

\vspace{0.25cm}

La tabella mostra che
\begin{enumerate}
\item $s_n$ non converge numericamente seppure lo faccia teoricamente, 
\item $t_n$ converge numericamente e teoricamente, qualora per calcolare $t_n$, $n < N$ si parta da un valore $m \gg N$. Nel nostro esempio, per $N=100$, siamo partiti da $m=500$.
\end{enumerate}

\vspace{0.25cm}

La motivazione di tale comportamento consiste nel fatto che la prima successione, cio\mye quella in avanti, amplifica enormemente l'errore in avanti mentre la seconda, detta all'indietro, fa esattamente il contrario.

In effetti anche se operassimo in aritmetica esatta, per quanto concerne la successione in avanti $I_n$, se $I^*_0=I_0+\epsilon$ allora la successione calcolata $I^*_{n+1}=1-(n+1)I^*_n$ \mye tale che  $I^*_n-I_n=n! \epsilon$, giustificando i risultati ottenuti.

\subsection{Esempio 4: sulla somma ${\frac {(1+\eta)-1} {\eta} }$}
Calcoliamo l'errore relativo tra $1$ e il valore che si ottiene
valutando in doppia precisione 
$$f(\eta)={\frac {(1+\eta)-1} {\eta} }$$   
con $\eta=10^{-1}, 10^{-2}, \ldots, 10^{-15}$. 
Ovviamente il risultato esatto di questa divisione \mye sempre $1$, per qualsiasi $\eta$. 


Notiamo che il calcolatore effettuer\mya invece
$f^*(\eta) $
ottenuta
\begin{enumerate}
\item valutando $a_1=fl(fl(1)+fl(\eta))$;
\item valutando $a_2=fl(fl(a_1)-fl(1))$;
\item valutando $a_3=fl(fl(a_2):fl(\eta))$.
\end{enumerate}
in cui ovviamente $fl(1)=1$, $fl(a_1)=a_1$, $fl(a_2)=a_2$ visto che $a_1$, $a_2$ sono numeri floating point.

I risultati ottenuti nel valutare $f^*(\eta) $ sono potenzialmente diversi da $f(\eta)=1$ (ma non per gli interi, perch\mye?).

Ricordando la notazione esponenziale in cui ad esempio {\tt{3.1e-02}}=$3.1 \cdot 10^{-2}$, e che l'errore relativo risulta
$$
\frac{|f(\eta)-f^*(\eta) |}{|f(\eta)|}=\frac{|1-f^*(\eta) |}{1}=|1-f^*(\eta)|
$$
abbiamo

\vspace{0.25cm}

\begin{center}
\begin{tabular}[h!]{| c | c | c |} 
\hline
$\eta$  & $f^*(\eta)$ & ${\mbox{rel.err.}}$ \\
\hline
 	 $1.00e-01$ & $1.000000000000001e+00$ & $8.881784e-16$ \\
 	 $1.00e-02$ & $1.000000000000001e+00$ & $8.881784e-16$ \\
 	 $1.00e-03$ & $9.999999999998899e-01$ & $1.101341e-13$ \\
 	 $1.00e-04$ & $9.999999999998899e-01$ & $1.101341e-13$ \\
 	 $1.00e-05$ & $1.000000000006551e+00$ & $6.551204e-12$ \\
 	 $1.00e-06$ & $9.999999999177334e-01$ & $8.226664e-11$ \\
 	 $1.00e-07$ & $1.000000000583867e+00$ & $5.838672e-10$ \\
 	 $1.00e-08$ & $9.999999939225290e-01$ & $6.077471e-09$ \\
 	 $1.00e-09$ & $1.000000082740371e+00$ & $8.274037e-08$ \\
 	 $1.00e-10$ & $1.000000082740371e+00$ & $8.274037e-08$ \\
 	 $1.00e-11$ & $1.000000082740371e+00$ & $8.274037e-08$ \\
 	 $1.00e-12$ & $1.000088900582341e+00$ & $8.890058e-05$ \\
 	 $1.00e-13$ & $9.992007221626409e-01$ & $7.992778e-04$ \\
 	 $1.00e-14$ & $9.992007221626409e-01$ & $7.992778e-04$ \\
 	 $1.00e-15$ & $1.110223024625157e+00$ & $1.102230e-01$ \\
\hline
%\caption{I valori assunti per valutare numericamente per certi $\eta$ la quantit\mya ${\frac {(1+\eta)-1} {\eta} }$.}
\end{tabular}
\end{center}

\vspace{0.25cm}

La perdita di precisione \mye dovuta essenzialmente al fenomeno di cancellazione, visto che $1+\eta \approx 1$. Ricordato che {\tt{eps}} \mye un numero macchina e che per numeri pi\myu piccoli $\eta$ in modulo di {\tt{eps}} si ha $1+x=1$ non sorprende inoltre sia

\vspace{0.25cm}

\begin{center}
\begin{tabular}[h!]{| c | c | c |} 
\hline
$\eta$  & $f^*(\eta)$ & ${\mbox{rel.err.}}$ \\
\hline
 ${\mbox{{\tt{eps}}}}$ & $1.000000000000000e+00$ & $0.000000e+00$ \\
 	 ${\mbox{{\tt{eps/2}}}}$ & $0.000000000000000e+00$ & $1.000000e+00$ \\
 	 ${\mbox{{\tt{eps/3}}}}$ & $0.000000000000000e+00$ & $1.000000e+00$ \\
\hline
%\caption{I valori assunti per valutare numericamente per certi $\eta$ la quantit\mya ${\frac {(1+\eta)-1} {\eta} }$.}
\end{tabular}
\end{center}

\vspace{0.25cm}

In effetti, per questi ultimi valori $\eta$ pi\myu piccoli in modulo di {\tt{eps}}, visto che comunque $fl(\eta) > 0$ (pensarci su!),
$$
{\frac {(1+fl(\eta))-1} {fl(\eta)} }={\frac {1-1} {fl(\eta)} }=0.
$$



 \subsection{Esempio 5: Sulla valutazione di $\tan(\arctan(x))=x$}
Siano $f:=\tan$ e $g:=\arctan$. Si consideri la funzione composta $$f 
\cdot 
g(x):=f(g(x))=\tan(\arctan(x))=x, \,\,x \in (-\infty,+\infty).$$ 

Ci domandiamo se l'implementazione numerica di questa funzione goda ancora di questa propriet\'a.

Valutiamo la funzione $g$ numericamente, mediante la corrispettiva versione numerica $g^*$ per
$$x=10^{k}, {\mbox{ per }}  k=-20+40h, {\mbox{ e }} 
h=0,0.1,0.2,\ldots,1$$ 
e valutiamo l'errore relativo compiuto, ovvero $|g(x)-g^*(x)|/|g(x)|$, notando che $g(x) \neq 0$ visto che nei casi in questione si ha $x \neq 0$.
Descriviamo questi risultati nella tabella che segue.

\vspace{0.25cm}

\begin{center}
\begin{tabular}[h!]{| c | c | c |} 
\hline
$x$  & $g^*(x)$ & ${\mbox{rel.err.}}$ \\
\hline
 	 $1.0e-20$ & $1.000000000000000e-20$ & $0.00000e+00$ \\
 	 $1.0e-16$ & $1.000000000000000e-16$ & $0.00000e+00$ \\
 	 $1.0e-12$ & $1.000000000000000e-12$ & $0.00000e+00$ \\
 	 $1.0e-08$ & $1.000000000000004e-08$ & $0.00000e+00$ \\
 	 $1.0e-04$ & $1.000000000000000e-04$ & $0.00000e+00$ \\
 	 $1.0e+00$ & $9.999999999999999e-01$ & $1.11022e-16$ \\
 	 $1.0e+04$ & $9.999999999990539e+03$ & $9.46056e-13$ \\
 	 $1.0e+08$ & $9.999999999542369e+07$ & $4.57631e-11$ \\
 	 $1.0e+12$ & $1.000071916854289e+12$ & $7.19117e-05$ \\
 	 $1.0e+16$ & $1.633123935319537e+16$ & $3.87677e-01$ \\
 	 $1.0e+20$ & $1.633123935319537e+16$ & $6.12223e+03$\\
 	 \hline
\end{tabular}
\end{center}

\vspace{0.25cm}

Osserviamo che i valori critici sono tra $1$ e $10^4$. In effetti un ulteriore analisi propone

\vspace{0.25cm}

\begin{center}
\begin{tabular}[h!]{| c | c | c |} 
\hline
$x$  & $g^*(x)$ & ${\mbox{rel.err.}}$ \\
\hline
$1.000000000000000e+00$ & $9.999999999999999e-01$ & $1.11022e-16$ \\
 	 $2.511886431509580e+00$ & $2.511886431509580e+00$ & $0.00000e+00$\\
 	 $6.309573444801933e+00$ & $6.309573444801929e+00$ & $7.03834e-16$\\
 	 $1.584893192461114e+01$ & $1.584893192461114e+01$ & $1.12081e-16$\\
 	 $3.981071705534973e+01$ & $3.981071705534988e+01$ & $3.56961e-15$\\
 	 $1.000000000000000e+02$ & $1.000000000000010e+02$ & $1.00897e-14$\\
 	 $2.511886431509580e+02$ & $2.511886431509614e+02$ & $1.38042e-14$\\
 	 $6.309573444801930e+02$ & $6.309573444802187e+02$ & $4.07210e-14$\\
 	 $1.584893192461114e+03$ & $1.584893192461186e+03$ & $4.54778e-14$\\
 	 $3.981071705534973e+03$ & $3.981071705533795e+03$ & $2.95849e-13$\\
 	 $1.000000000000000e+04$ & $9.999999999990539e+03$ & $9.46056e-13$\\
 	 \hline
\end{tabular}
\end{center}

\vspace{0.25cm}

Quindi numericamente $g$ e $g^*$ non coincidono, anche per valori non eccessivi di $x$.


















\section{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 {\sl{bit}} una cifra nel sistema (usualmente binario) usato {\cite[p. 48]{COM}}, {\cite[p. 47]{MATNUM}}. 

Nello standard IEEE 754, i numeri in virgola mobile a  precisione singola  sono memorizzati in un vettore
detto {\sl{parola}} di 32 bits, mentre quelli a precisione doppia occupano 
due parole consecutive da 32 bits. Un numero {\sl{non nullo-normalizzato}} si scrive come 
\begin{equation}\label{1aa}
x = (-1)^S \cdot 2^{E-Bias} \cdot (1+F)
\end{equation}
dove $S=0$ oppure $S=1$ determina il {\sl{segno}}, $E-Bias$ \mye 
l'{\sl{esponente}} cui va sottratto il valore del BIAS ed $F$ la 
{\sl{mantissa}} con $F=\sum^{+\infty}_{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 {\sl{numeri macchina}}

\begin{equation}\label{1}
x = (-1)^S \cdot 2^{E-Bias} \cdot (1+F)
\end{equation}
con $E$ intero tale che $|E| < 2^{N_E}$ per un certo prefissato ${{N_E}}$ e  $F=\sum^{N_F}_{k=1}a_k 2^{-k}$,



Alcuni esempi
{\smsp}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
Precisione & $N_E$ & $N_F$ & {{Bias}} \\ 
\hline
Singola    & 8           & 23      & 127  \\
Doppia     & 11        & 52       & 1023   \\
\hline
\end{tabular}
\end{center}
{\smsp}
\noindent


Partiamo con un esempio semplice. Il numero $12$ corrisponde in singola precisione a
$$[0 \,|\, 10000010 \,|\, 10000000000000000000000]$$

\noindent
Cerchiamo di comprenderne il significato:
{\smsp}
\begin{enumerate}
\item $0$ determina il segno. Poich\mye $(-1)^0=1$ il segno \mye $+$.
{\smsp}
\item $10000010$ determina il valore della potenza di $2$. Essendo il BIAS uguale a 127 e (si 
legga da destra a sinistra!) 
$E=10000010=0 \cdot 2^0+1 \cdot 2^1+0 \cdot 2^2+0 \cdot 2^3+0 \cdot 2^4+0 \cdot 2^5+0 \cdot 2^6+1 \cdot 2^7=2+128=130$, 
l'esponente vale $130-127=3$. Quindi $2^{E-Bias}=2^3=8$.
{\smsp}
\item $F=10000000000000000000000$ \mye la componente relativa alla 
mantissa. La quantit\mya, letta da sinistra a destra, riguarda la parte 
decimale $F$ in (\ref{1}). Quindi $F=1 \cdot 2^{-1}+0 \cdot 2^{-2}+\ldots=0.5$. Ora 
($1+F$) si legge come $1+F=1.5$.
\end{enumerate}

\vspace{0.5cm}
In altre parole la rappresentazione di $12$ \mye esatta ed il numero \mye 
rappresentato come $12=8 \cdot 1.5$.

Non convinti, proviamo un altro esempio, ad esempio 126 che corrisponde in 
singola precisione a 
$$[0 \,|\, 10000101\,|\, 11111000000000000000000]$$
{\smsp}
\begin{enumerate}
\item $0$ determina il segno. Poich\mye $(-1)^0=1$ il segno \mye $+$.
{\smsp}
\item $10000101$ determina il valore. Essendo il BIAS uguale a 127 e (si 
legga da destra a sinistra!) 
$E=10000101=1 \cdot 2^0+0 \cdot 2^1+1 \cdot 2^2+0 \cdot 2^3+0 \cdot 2^4+0 \cdot 2^5+0 \cdot 2^6+1 \cdot 2^7=1+4+128=133$, 
l'esponente vale $133-127=6$. Quindi $2^{E-Bias}=2^6=64$. Viene da chiedersi perch\mye introdurre il BIAS. La risposta \mye che rappresenta un alternativa all'introduzione di un bit di segno per l'esponente. 
{\smsp}
\item $F=11111000000000000000000$ \mye la componente relativa alla 
mantissa. La quantit\mya, letta da sinistra a destra, riguarda la parte 
decimale $F$ in (\ref{1}). Quindi $F= 1 \cdot 2^{-1}+1 \cdot 2^{-2}+ 
1 \cdot 2^{-3}+1 \cdot 2^{-4}+1 \cdot 2^{-1}+\ldots  
=0.5+0.25+0.125+0.0625+0.03125=0.96875$. Ora ($1+F$) si legge come 
$1+F=1.96875$.
\end{enumerate}

\vspace{0.5cm}
In altre parole la rappresentazione di $126$ \mye esatta ed il numero \mye 
rappresentato come $126=64\,\cdot\,1.96875$.

La distanza tra un numero 
in virgola mobile $x^-$ ed il successivo $x^+$, \mye 
$$2^{-N_F}  \,\cdot\,  2^{{E-bias}}.$$
Per convincersene, supponiamo che siano 
$$x^-= (-1)^S \cdot 2^{E-Bias} \cdot (1+F^-), \,\, x^+= (-1)^S \cdot 2^{E-Bias} \cdot (1+F^+).$$
Allora, essendo $F^+-F-=2^{-N_F}$ (pensarci su!) facilmente $$|x^+-x^-|=2^{E-Bias} \cdot (F^+-F-)=2^{-N_F}  \,\cdot\,  2^{{E-bias}}.$$

La cardinalit\mya dei numeri macchina in IEEE 754 \mye piuttosto elevata. Ad esempio, in dop\-pia pre\-ci\-sio\-ne sono $$2^{64} \approx 1.844674407370955e+019.$$
La {\sl{precisione di macchina}}, cio\mye 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 ${\tt{eps}}$ in Matlab 7.6.0 vale

\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}
























\subsection{Facoltativo: Su {\tt{eps}}, {\tt{realmax}}, {\tt{realmin}}, overflow e underflow in IEEE754}
I numeri nella rappresentazione floating-point stanno in un intervallo 
limitato, quindi l'{\sl{infinitamente piccolo}} e l'{\sl{infinitamente 
grande}} non sono rappresentabili come numeri in virgola mobile, nel 
calcolatore. Quando si raggiunge un valore talmente piccolo da non essere 
pi\myu distinto dallo zero, si parla di {\sl{underflow}}, mentre quando si 
eccede il massimo numero rappresentabile, si parla di {\sl{overflow}}.
Il limite inferiore (che non sia {\sl{denormalizzato}}!) \mye pari a $2^{ 1-{{bias}} } $ ed in Matlab \mye 
espresso dalla variabile {\tt{realmin}}.
Il limite superiore \mye pari a $2^{2^N_E-2-bias} \cdot (1+(1- 
2^{\mbox{-Fbits}} ))$ ed in Matlab \mye espresso dalla variabile 
{\tt{realmax}}.
Verifichiamolo.
{\smsp}
\begin{enumerate}
\item Dalla tabella, in precisione doppia, si ha che {\sl{bias}}=1023 e 
quindi $$2^{ 1-1023} =2^{ -1022 }=2.225073858507201e-308 $$
Se digitiamo nella shell di Matlab il comando {\tt{realmin}} otteniamo 
proprio questo valore.
{\smsp}
\item Dalla tabella, in precisione doppia, si ha che {\sl{bias}}$=1023$, 
$N_E=11$, ed $N_F=52$. Conseguentemente
\begin{eqnarray}
{\mbox{realmax}} & = & 2^{2^{N_E}-2-bias} \, \cdot \, (1+(1- 2^{N_F} )) \nonumber\\  
& = & 2^{2^{11}-2-1023} \, \cdot \, (1+(1- 2^{-52} )) \nonumber \\
& = & 2^{2^{11}-2-1023} \, \cdot \, (2- 2^{-52}) = \nonumber\\
& = & 1.797693134862316e+308 \nonumber
\end{eqnarray}

\noindent
Per convincercene, digitiamo nella shell di Matlab/Octave
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}
Quindi digitando {\tt{realmax}} nella shell di Matlab otteniamo proprio questo 
valore.
\end{enumerate}

Osserviamo che in realt\mya {\tt{realmin}} non sembra essere il numero pi\myu piccolo in valore assoluto. Testiamolo sulla shell di Matlab.
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}
\noindent
Siccome in un costrutto logico il valore $1$ corrisponde a {\tt{true}}, si ricava che 
$${\mbox{realmin}}/(2^{52}) \approx 4.940656458412465{\mbox{e}}-324$$ \mye ancora un numero macchina, ma che fa parte dei cosidetti {\sl{numeri denormalizzati}} (cf. {\cite{MATNUM}}, p.49) e cio\mye del tipo 
\begin{equation}\label{1a}
x = (-1)^S \, \cdot \, 2^{E-Bias} \, \cdot \, (0.F).
\end{equation}

Per curiosit\mya vediamo il codice Matlab che implementa {\tt{realmin}}:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}
Per quanto concerne {\tt{pow2}}:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

>> 

\end{Verbatim}
























\section{Alcuni esercizi in Matlab}

\vspace{0.5cm}


\begin{enumerate}

% ESERCIZIO 1


\item {\bf{Esercizio facoltativo difficile}}. Esistono pi\myu successioni che approssimano $\pi$ e sono attribuite ad Archimede. Ad esempio, osservato che questi non \mye altro che l'area del cerchio avente raggio $1$, si inscrivono nel cerchio al variare di $p=2,3,\ldots$ dei poligoni regolari aventi $2^p$ lati. Osserviamo che per $p=2$, si deve determinare l'area del quadrato inscritto, che con facili conti \mye 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).

 \begin{figure}\label{fig1s}
 \centering
   {\includegraphics[scale=0.35,clip]{archimede_8.eps}}
   {\includegraphics[scale=0.35,clip]{archimede_16.eps}}
 \caption{Grafico che illustra le suddivisioni considerate dall'algoritmo di Archimede, per $p=3$, $p=4$.}
 \end{figure}

Se indichiamo con $\theta$ l'angolo al centro di ogni triangolo, si verifica che la sua area \mye $\frac{\sin {(\theta)}}{2}$. In particolare per un poligono di $2^p$ lati, si ha $\theta_p=\frac{2\pi}{2^p}$ e l'area del poligono inscritto, costituito da $2^p$ triangoli uguali, \mye 
$$I_p=2^p \, \frac{\sin {(\theta_p)}}{2} = 2^p \, \frac{\sin {(\frac{2\pi}{2^p})}}{2}.  $$
In realt\mya Archimede us\myo questa ed altre idee per cercare di approssimare $\pi$ con stime dall'alto e dal basso, via anche poligoni circoscritti. Come citato in {\cite{WIKAR}}:

\vspace{0.25cm}

\begin{quote}
Nel breve lavoro "{\sl{La misura del cerchio}}'' viene dimostrato anzitutto che un cerchio \mye equivalente a un triangolo con base eguale alla circonferenza e altezza eguale al raggio. 

Tale risultato \mye 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\myo approssimare arbitrariamente il rapporto tra circonferenza e diametro di un cerchio dato, rapporto che oggi si indica con $\pi$. Le stime esplicitamente ottenute limitano questo valore fra 22/7 e 3 + 10/71..
\end{quote}

\vspace{0.25cm}

Secondo {\cite{WOL}}, Archimede us\myo un altro algoritmo per approssimare $\pi$. 

Visto che $2\pi $ non \mye altro che la lunghezza della circonferenza del cerchio avente raggio uguale a $1$ e questa \mye maggiore del perimetro $a_k$ di un qualsiasi poligono regolare di $n=6 \cdot 2^k$ lati inscritto e minore del perimetro di un qualsiasi poligono regolare di $n=6 \cdot 2^k$ lati circoscritto $b_k$, misurando $a_k$ e $b_k$ si pu\myo affermare che $$a_k \leq 2\pi \leq b_k. $$

\noindent
Si pu\myo inoltre provare che 
$$a_k=2n \,  \sin {\left(\frac{\pi}{n} \right)},$$
$$b_k=2n \, \tan {\left(\frac{\pi}{n}\right)}.$$
Stimare $\pi$ come $$a_k \leq 2 \, \pi \leq b_k.$$ 

\noindent
Per quale $n$ si riesce a determinare $\pi$ con $10$ cifre decimali esatte? Pensare se sia giusto dire $b_k-a_k < 10^{-10}$ allora $\pi$ \mye approssimato da $a_k$ con 10 cifre decimali esatte.


\noindent
{\bf{A.}} Si approssimi $\pi$ 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.

\noindent
{\bf{B.}} Stimare $\pi$ come $$a_k \leq \pi \leq 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 $\pi$ con $10$ cifre decimali esatte? Per $k=0,\ldots,4$ si ha
\begin{eqnarray}
3.00000 & \leq \pi \leq & 3.46410 \\	
3.10583 & \leq \pi \leq & 3.21539 \\
3.13263 & \leq \pi \leq & 3.15966 	\\
3.13935 & \leq \pi \leq & 3.14609 	\\
3.14103 & \leq \pi \leq & 3.14271
\end{eqnarray}

\begin{remark} 
Il fatto che ci sia un $\pi$ nei secondi membri di alcune espressioni non \mye {\sl{un gatto che si morde la coda}}. Infatti \mye solo una trascrizione in radianti di un angolo che ha un significato geometrico indipendente dalla quantit\mya $\pi$. In questi esercizi, per semplificare la questione si usa per\myo $\pi$ per calcolare $\pi$, cosa che dimostra una volta in pi\myu l'abilit\mya di Archimede, che calcol\myo tali quantit\mya con metodi elementari e geometrici.
\end{remark} 

\vspace{0.25cm}

Si cita il fatto che per $k=4$ Archimede riusci' ad affermare che 
$$\frac{223}{71} \leq \pi \leq \frac{22}{7}.$$
{\smsp}
% ESERCIZIO 5.
\item {\bf{Esercizio, facile facoltativo}}. Posto $h=10^{-1},10^{-2},\ldots,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$?
{\smsp}
% ESERCIZIO 6.
\item {\bf{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.
{\smsp}
% ESERCIZIO 7.
\item {\bf{Esercizio, facile facoltativo}}. Fissato $\theta=0:(2\pi / 1001):2 \pi$, si valuti $\sin^2{(\theta)} + \cos^2{(\theta)}$ e si calcoli l'errore assoluto rispetto il valore $1$.
\end{enumerate}



























\section{Online}
Alcuni siti utili per la comprensione della lezione:
{\smsp}
\begin{enumerate}
\item \href{http://it.wikipedia.org/wiki/Archimede}{http://it.wikipedia.org/wiki/Archimede}
{\smsp}
\item \href{http://utenti.quipo.it/base5/numeri/pigreco.htm}{http://utenti.quipo.it/base5/numeri/pigreco.htm}
{\smsp}
\item \href{http://it.wikipedia.org/wiki/Pi\_greco}{http://it.wikipedia.org/wiki/Pi\_greco}
{\smsp}
\item \href{http://it.wikipedia.org/wiki/Teorema\_fondamentale\_del\_calcolo\_integrale}{http://it.wikipedia.org/wiki/Teorema\_fondamentale\_del\_calcolo\_integrale}
{\smsp}
\item \href{http://it.wikipedia.org/wiki/Floating\_point}{http://it.wikipedia.org/wiki/Floating\_point}
{\smsp}
\item \href{http://it.wikipedia.org/wiki/IEEE\_754r}{http://it.wikipedia.org/wiki/IEEE\_754r}
{\smsp}
\item \href{http://it.wikipedia.org/wiki/IEEE\_754}{http://it.wikipedia.org/wiki/IEEE\_754}
{\smsp}
\item \href{https://en.wikipedia.org/wiki/IEEE\_754-1985}{https://en.wikipedia.org/wiki/IEEE\_754-1985}
{\smsp}
\item \href{http://mathworld.wolfram.com/ArchimedesAlgorithm.html}{http://mathworld.wolfram.com/ArchimedesAlgorithm.html}
\end{enumerate}

\begin{thebibliography}{99} 
\bibitem{ATK} K. Atkinson, {\em An Introduction to Numerical Analysis\/}, Wiley, (1989).
\bibitem{ATKHAN} K. Atkinson and W. Han, {\em Elementary Numerical Analysis\/}, third edition, Wiley, (2004).
\bibitem{BOth} N. Brisebarre and collaborators {\em {Floating-Point Arithmetic}}, 2009.\\ {\href{https://perso.ens-lyon.fr/jean-michel.muller/chapitre1.pdf}{https://perso.ens-lyon.fr/jean-michel.muller/chapitre1.pdf}}
\bibitem{COM} V. Comincioli, {\sl{Analisi Numerica, metodi modelli applicazioni}}, Mc Graw-Hill, 1990. 
\bibitem{COMESE} V. Comincioli, {\sl{Problemi di analisi numerica}}, Mc Graw-Hill, 1991.
\bibitem{FRO} M. Frontini e E. Sormani, {\sl{Fondamenti di Calcolo Numerico, problemi in laboratorio}}, Apogeo, 2005.
\bibitem{GOLD} D. Goldberg, {\sl{What Every Computer Scientist Should Know About Floating-Point Arithmetic}},  Association for Computing Machinery, Inc, Computing Surveys, March 1991.
\bibitem{HUC} T. Huckle , {\href{http://wwwzenger.informatik.tu-muenchen.de/persons/huckle/bugse.html}{http://wwwzenger.informatik.tu-muenchen.de/persons/huckle/bugse.html}}. 
\bibitem{MF} J.H. Mathews e K.D. Fink, {\sl{Numerical Methods using Matlab}}, Prentice Hall, 1999. 
\bibitem{OBI} Network Theory, {\sl{GNU Octave Manual Version}}, {\href{http://www.network-theory.co.uk/docs/octave3/octave\_160.html}{http://www.network-theory.co.uk/docs/octave3/octave\_160.html}}. 
\bibitem{QUAold} A. Quarteroni, {\sl{Elementi di calcolo numerico}}, seconda edizione, Progetto Leonardo, 1999. 
\bibitem{QUA} A. Quarteroni e F. Saleri, {\sl{Introduzione al calcolo scientifico}}, Springer Verlag, 2006. 
\bibitem{MATNUM} A. Quarteroni, R. Sacco e F. Saleri, {\sl{Matematica Numerica}}, Springer Verlag, 1998. 
\bibitem{STO} J. Stoer, {\sl{Introduzione all'analisi numerica}}, Zanichelli, 1984. 
\bibitem{MOL} The MathWorks Inc., {\sl{Numerical Computing with Matlab}},{\href{http://www.mathworks.com/moler}{http://www.mathworks.com/moler}}.
\bibitem{TRENUM} 
L.N. Trefethen, {\sl{Numerical Analysis}},{\href{http://people.maths.ox.ac.uk/trefethen/NAessay.pdf}{http://people.maths.ox.ac.uk/trefethen/NAessay.pdf}}.
\bibitem{UNIPO} Web Page , \\ {\href{http://utenti.quipo.it/base5/numeri/pigreco.htm}{http://utenti.quipo.it/base5/numeri/pigreco.htm}}.
\bibitem{WIKAR} Wikipedia (Archimede), \\ {\href{http://it.wikipedia.org/wiki/Archimede}{http://it.wikipedia.org/wiki/Archimede}}
\bibitem{WIKCI} Wikipedia (Teorema Fondamentale del Calcolo Integrale), \\ \href{http://it.wikipedia.org/wiki/Teorema\_fondamentale\_del\_calcolo\_integrale}{http://it.wikipedia.org/wiki/Teorema\_fondamentale\_del\_calcolo\_integrale}
\bibitem{WIKFL} Wikipedia (Floating Point), \\ \href{http://it.wikipedia.org/wiki/Floating\_point}{http://it.wikipedia.org/wiki/Floating\_point}
\bibitem{WIKFL2} Wikipedia (IEEE 754), \\ \href{http://it.wikipedia.org/wiki/IEEE\_754}{http://it.wikipedia.org/wiki/IEEE\_754}
\bibitem{WIKFL3} Wikipedia (IEEE 754r), \\ \href{http://it.wikipedia.org/wiki/IEEE\_754r}{http://it.wikipedia.org/wiki/IEEE\_754r}
\bibitem{WIKPI} Wikipedia (Pi Greco), \\ {\href{http://it.wikipedia.org/wiki/Pi\_greco}{http://it.wikipedia.org/wiki/Pi\_greco}}
\bibitem{WOL} Wolfram (Archimedes' Algorithm), \\ \href{http://mathworld.wolfram.com/ArchimedesAlgorithm.html}{http://mathworld.wolfram.com/ArchimedesAlgorithm.html}
\end{thebibliography}
\end{document}





\end{document}



\section{Alcune frasi celebri}

In {\cite{UNIPO}}, oltre a varie curiosit\mya su $\pi$ si cita un curioso aneddoto riguardante il noto matematico J. H. Conway:


``{\sl{Un giorno decisi di imparare a memoria le prime mille cifre del pi greco}} - ricorda Conway - {\sl{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\mya 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}}''.


``{\sl{E' stato divertente}} - continua Conway - {\sl{perch\mye 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\myi di seguito, alternandoci nella recita: in questo modo siamo arrivati a memorizzare le mille cifre del pi greco}}''.


Pure la morte di Archimede \mye coperta dalla leggenda {\cite{WIKAR}}:

{\sl{Ad un tratto entr\myo nella stanza un soldato e gli ordin\myo 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\myo, sguain\myo la spada e lo uccise}}.

\noindent
Frasi celebri attribuite ad Archimede 
{\smsp}
\begin{itemize}
\item {\sl{Noli turbare circulos meos}} cio\mye {\sl{Non turbare i miei cerchi}}.
{\smsp}
\item le sue ultime parole pare furono {\sl{Soldato, stia lontano dal mio disegno}}.
\end{itemize}




\subsection{Implementazione in Matlab di una certa equazione di secondo grado}
Dato il polinomio di secondo grado $x^2 + 2\,p x -q$, con 
${\sqrt{p^2+q}}\geq 0$ calcoliamo la 
radice 
\begin{equation}\label{quadeq}
y=-p+{\sqrt{p^2+q}}.
\end{equation}
mediante due algoritmi (uno stabile e uno instabile)
\begin{equation}\label{quadeq2R}
y=-p \pm {\sqrt{p^2+q}}.
\end{equation}
e
\begin{equation}\label{quadeqrR}
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}})} }.
\end{equation}

Scriviamo un programma {\tt{radicesecgrado.m}} in Matlab che illustri i due algoritmi.
{\vspace{0.25cm}}
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

% 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

\end{Verbatim} 
\noindent
Digitiamo quindi da shell Matlab/Octave il comando {\tt{radicesecgrado}} e otteniamo 
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}
\noindent
Come previsto, il secondo algoritmo si comporta notevolmente meglio del primo, che compie un errore relativo dell'ordine di circa $10^{-9}$. 

{\vspace{0.25cm}}

\noindent
{\bf{Esercizio veloce}}. Testare il codice {\tt{radicesecgrado.m}} per l'esempio in cui
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}


\subsection{Esempio 2: Approssimazione numerica di $\pi$}

Storicamente sono state scoperte diverse successioni convergenti sempre pi\myu rapidamente a $\pi$ (cf. {\cite{WIKPI}}).
In questa sezione ci interesseremo a 3 di queste mostrando sia come le velocit\mya di convergenza possano essere diverse, sia come possano insorgere questioni di instabilit\mya {\cite[p. 16]{FRO}}.

\noindent
Si implementino le successioni $\{u_n\}$, $\{z_n\}$, definite 
rispettivamente come
$$\left\{ \begin{array} {l}  
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{array} \right.$$
e
$$\left\{ \begin{array} {l}  
z_1=1, \, \, z_2=2 \\
z_{n+1} = 2^{n- { \frac {1} {2} } }      { \sqrt { 1 - \sqrt{1-4^{1-n} 
\cdot z_n^2}   } } 
\end{array} \right.$$
che {\sl{teoricamente}} convergono a $\pi$.
Si implementi poi la successione, diciamo $\{y_{n}\}$, che si ottiene 
{\sl{razionalizzando}}, cio\mye moltiplicando numeratore e denominatore 
per 
$$ { \sqrt { 1 + \sqrt{1-4^{1-n} \cdot z_n^2}   } } $$
e si calcolino $u_m$, $z_m$ e $y_m$ per $m=2,3,\ldots,40$ (che teoricamente dovrebbero approssimare $\pi$).

\noindent
Si disegni in un unico grafico l'andamento dell'errore relativo di $u_n$, 
$z_n$ e $y_n$ rispetto a $\pi$. A tal proposito 
ci si aiuti con l'help di Matlab relativo al comando {\tt{semilogy}}. I 
grafici devono avere colori o patterns 
diversi.

{\vspace{0.25cm}}

\noindent
{\bf{Facoltativo}}: in un riquadro mettere un legame tra colore (e/o pattern) e 
successione, usando il comando Matlab {\tt{legend}} (aiutarsi con l'help). 
Ricordiamo che tale comando non esiste in vecchie versioni di GNU Octave. Un esempio del suo utilizzo in Octave \mye (cf. {\cite{OBI}}
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

legend ("sin (x)");

\end{Verbatim}
{\vspace{0.25cm}}

\subsection{Risoluzione} 

In seguito scriviamo un'implementazione di quanto richiesto commentando i 
risultati. Si salvi in un file {\tt{pigreco.m}} il codice
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

% 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.
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;

\end{Verbatim}
 \begin{figure}\label{fig1}
 \centering
   {\includegraphics[scale=0.5,clip]{plot_pigreco.eps}}
 \caption{Grafico che illustra le 3 successioni, rappresentate rispettivamente da {\tt{.}}, {\tt{+}} e {\tt{o}}.}
 \end{figure}


\subsubsection{Commenti alla Risoluzione}
{\smsp}
\begin{enumerate}
\item Non tutti i programmi sono functions. Alcuni sono semplicemente un 
codice che viene interpretato da Matlab (il cosidetto {\sl{programma 
principale}}). Usiamo funzioni solo quando vogliamo introdurre porzioni di 
codice che possono tornare utili a pi\myu programmi principali, o che 
semplificano la loro lettura.
{\smsp}
\item Pi\'u assegnazioni possono essere scritte in una stessa riga. Per convincersene si osservi la prima riga del file {\tt{pigreco.m}} dopo i 
commenti.
{\smsp}
\item L'istruzione descritta dopo {\tt{for n=2:40}} non richiede l'incremento della 
variabile $n$, in quanto ci\myo \mye automatico.
{\smsp}
\item Non serve il {\tt{;}} dopo l'{\tt{end}} che chiude il ciclo 
{\tt{for}}.   
{\smsp}
\item Con un po' di tecnica il ciclo {\tt{for}} \mye sostituibile col un 
ciclo {\tt{while}}. In questo caso per\myo bisogna incrementare la variabile $n$.
{\smsp}
\item Nel denominatore riguardante l'errore relativo scriviamo {\tt{pi}} e 
non {\tt{abs(pi)}} in quanto $\pi=|\pi|$. 
{\smsp}
\item Nel comando {\tt{fprintf}}, utilizziamo {\tt{$\backslash$n}} che manda a capo e {\tt{$\backslash$t}} che 
esegue un 
{\sl{tab}} 
(uno spazietto in avanti). Se non si scrive {\tt{$\backslash$n}}, Matlab 
scriver\mya di seguito l'output sulla stessa riga, rendendo difficile la lettura dell'errore 
relativo.
{\smsp}
\item Nel comando {\tt{fprintf}} alcune parti vengono scritte come testo 
altre come contenuto di una variabile. Consideriamo ad esempio gli
{\tt{fprintf}} nella prima successione:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}
Dopo essere andati a capo e spaziato a destra, Matlab scrive sul monitor 
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]
 
{\tt {[SEQ.1][INDEX]: } } 

\end{Verbatim}
e quindi valuta {\tt {n+1}} che viene stampato 
con {\sl{tre cifre decimali prima della virgola}} in notazione decimale. 
Scrive poi sul monitor {\tt {[REL.ERR.]:}}, e accede alla {\sl{cella di 
memoria}}
della variabile {\sl{relerru}} di cui considera la componente $n+1$-sima.
Di seguito stampa su monitor il suo valore con {\sl{due cifre decimali 
prima della 
virgola}}, {\sl{due cifre decimali dopo la virgola}} in notazione 
esponenziale.  
{\smsp}
\item Il comando {\tt{semilogy}} ha come primo argomento l'ascissa (che in 
questo caso sono gli indici di iterazione {\tt {1:41}}) e quale secondo 
argomento l'ordinata {\tt {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$, $\ldots$, 
dim( {\tt{relerru}}). Il grafico viene {\sl {tenuto}} grazie al comando 
di {\tt {hold on}} e di seguito si ripete il procedimento per {\tt {relerrz}} e {\tt {relerry}}. Si osservi 
che i vettori {\sl {ascissa}} e {\sl {ordinata}} devono essere della 
stessa dimensione e dello stesso tipo (cio\mye entrambi vettori riga o 
entrambi vettori colonna).
{\smsp}
\item Dall'help di {{semilogy}} si evince che se la variabile {\sl {ascissa}} non viene 
scritta allora Matlab indicizza automaticamente col vettore di naturali da 
1 alla dimensione del vettore {\sl {ordinata}}.
\end{enumerate}

\noindent
Per il risultato del plot si consideri la prima figura. Abbiamo indicato la prima successione con {\tt{.}}, la 
seconda con {\tt{+}} e la terza successione con {\tt{o}}. Osserviamo che in Octave, i grafici vengono rappresentati differentemente in quanto invece di un cerchietto viene talvolta disegnato un rombo.

\noindent
Dal punto di vista dell'analisi numerica si vede che 
{\smsp}
\begin{enumerate}
\item La prima successione converge molto lentamente a $\pi$, la seconda 
diverge mentre la terza converge velocemente a $\pi$. 
{\smsp}
\item Per alcuni valori $\{z_n\}$ e $\{y_n\}$ coincidono per alcune 
iterazioni per poi 
rispettivamente divergere e convergere a $\pi$. Tutto ci\myo \mye naturale 
poich\mye le 
due sequenze sono analiticamente (ma non numericamente) equivalenti.
{\smsp}
\item 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).
\end{enumerate}

\subsection{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 {\cite[Esercizio 1.9, p.35]{QUA}}. Problemi simili con una precisa discussione della propagazione dell'errore sono trattati pure in {\cite[p. 23, problema 11]{FRO}}

Sia 
\begin{equation}\label{321}
I_n=e^{-1} \int_0^1 x^n \, e^x \, dx
\end{equation}

\noindent
Per $n=0$ si ha $$I_0=e^{-1} \int_0^1 \, e^x \, dx=e^{-1}(e^{1}-1).$$
Per $n=1$, \mye facile verificare che, essendo
$$\int x \, \exp(x) dx=(x-1)\, \exp(x)+c,$$
dal secondo teorema del calcolo integrale (cf. {\cite{WIKCI}})
si ha che $I_1=e^{-1}$ e pi\myu in generale integrando per parti che
\begin{equation}\label{322}
I_{n+1}=e^{-1} \left( x^{n+1} \, e^x \mid_0^1 - (n+1) \int_0^1 x^n \, e^x 
\, dx  \right) =1-(n+1)\, I_n. 
\end{equation}

\noindent
Da (\ref{321}) essendo l'integranda $x^n \exp{x} >0 $ per $x \in (0,1]$
si ha $$I_n >0.$$ Inoltre $x \leq 1$ implica $x^2 \leq x$ e pi\myu in generale $x^{n+1} \leq x^n$ da cui $x^{n+1} \exp(x) \leq x^n \exp(x)$ e quindi $$I_{n+1} \leq I_n.$$
La successione $I_n$ \mye positiva e non crescente e quindi ammette limite $L$ finito.
Da ({\ref{322}}), calcolando il limite $L$ per $n \rightarrow \infty$ 
ad ambo i membri si ottiene portando $1$ a primo membro
$$L-1=-\lim_n  (n+1)\, I_n.$$
L'unica possibilit\mya affinch\mye $\lim_n  (n+1)\, I_n$ esista e sia finito \mye 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=+\infty,$$ il che \mye assurdo essendo $L < + \infty$ e quindi $L-1 < +\infty$.
{\smsp}
\begin{enumerate}
\item Si calcoli $I_n$ per $n=1,\ldots, 99$ mediante la successione {\sl{in avanti}}
$$\left\{ \begin{array}{l}  
s_1=e^{-1}\\
s_{n+1} = 1-(n+1) \, s_{n}
\end{array} \right.$$
{\smsp}
\item Fissato $m=500$, si calcoli la successione {\sl{all'indietro}}
$\{t_n\}_{n=1,\ldots,100}$ definita come 
$$\left\{ \begin{array}{l}  
t_{2m}=0\\
t_{n-1} = {\frac {1-t_n} {n}}
\end{array} \right.$$
Si osservi che per raggiungere tale obiettivo bisogna calcolare i termini 
$$t_m, \,t_{m-1},\ldots,t_{100},t_{99},\ldots,t_2,t_1.$$
{\bf{Nota}}: la successione all'indietro deriva dall'osservare che per $n$ sufficientemente grande $I_n \approx 0$. Quindi, posto per $n$ sufficientemente grande $I_n=0$ riscrivendo la successione in avanti come successione all'indietro (cio\mye $I_n$ in funzione di $I_{n+1}$), abbiamo 
$$I_{n} = {\frac {1-I_{n+1}} {n+1}}.$$
{\smsp}
\item Si disegni in un unico grafico semi-logaritmico (usare 
{\tt{semilogy}} e 
{\tt{subplot}}) l'andamento di $|s_n|$ e  
$|t_n|$ per $n=1,\ldots,100$.
{\smsp}
\item 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$.
{\smsp}
\item (Esercizio non banale). Si calcolino a mano 
$e^{(s)}_m:=I_m-s_m$ in funzione di $e^{(s)}_1$ e $e^{(t)}_1:=I_1-t_1$ in 
funzione di 
$e^{(t)}_{m}$. Infine si spieghi l'andamento oscillante della successione 
$\{s_n\}$.

\noindent
Suggerimento: si scriva $e^{(s)}_1=\delta$ (o equivalentemente $s_1=e^{-1}+\delta$) ed in seguito si esprima $e^{(s)}_m:=I_m-s_m$ in funzione di $\delta$.
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$. 

\end{enumerate}

{\vspace{0.25cm}}

\subsubsection{Risoluzione}
Di seguito vediamo un'implementazione di quanto richiesto. L'ultimo punto del problema lo lasciamo al lettore.
Scriviamo il codice in un file {\tt{succricorrente.m}}:

{\vspace{0.25cm}}
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim} 
\subsubsection{Commenti alla Risoluzione}

 \begin{figure}\label{fig2}
 \centering
   {\includegraphics[scale=0.5,clip]{ricorrente.eps}}
 \caption{Grafico che illustra i valori assoluti assunti dalla successione in avanti (in nero) e all'indietro (in rosa magenta)}
 \end{figure}

\noindent
Alcune osservazioni sul codice
{\smsp}
\begin{itemize}
\item abbiamo usato un comando del tipo {\tt {for n=M:-1:2}} cosicch\mye il ciclo for parte da
$M$ e arriva a 2, sottraendo di volta in volta $1$;
\item abbiamo inizializzato con un comando del tipo {\tt {zeros(M,1)}} il
vettore colonna $t$, che altrimenti non \mye in memoria (e quindi le componenti pi\myu piccole di
{\tt{t(n-1)}} sarebbero inaccessibili poich\mye indefinite);
{\smsp}
\item abbiamo applicato un grafico semilogaritmico poich\'e quello usuale non dice granch\mye (sperimentarlo).
{\smsp}
\item in figura abbiamo plottato le successioni $\{|s_i|\}$, $\{|t_i|\}$ e non $\{s_i\}$, $\{t_i\}$. Questo non \mye un problema in quanto 
$$s_i \rightarrow 0 \Leftrightarrow |s_i| \rightarrow 0$$
$$t_i \rightarrow 0 \Leftrightarrow |t_i| \rightarrow 0.$$
D'altro canto la prima successione diverge assumendo anche valori negativi, rendendo difficile la comprensione del grafico.
{\smsp}
\item Il comando {\tt{subplot}} permette di plottare pi\myu grafici nella stessa finestra di Matlab/Octave, senza sovrapporli. Per ulteriori informazioni, dall'help di Matlab:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

 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.

\end{Verbatim}
\end{itemize}

\noindent
Numericamente, dopo aver digitato sulla shell di Matlab/Octave {\tt{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, \mye evidente che la prima non fornisce la soluzione, mentre la seconda sembra convergere a $0$.

\noindent
Il secondo grafico (si veda la Figura 5) per vari $M=2 \cdot 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\myu accurata. In particolare l'errore assoluto compiuto \mye stampato 
nel display dal programma {\tt{succricorrente.m}} come segue:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

>> 

\end{Verbatim} 
 \begin{figure}\label{fig3}
 \centering
   {\includegraphics[scale=0.5,clip]{ricorrente2.eps}}
 \caption{Grafico che illustra il valore assoluto della successione
all'indietro nell'approssimare $\exp{(-1)}$.}
 \end{figure}

\section{Complessit\mya computazionale}
Ricordato che un {\sl{algoritmo}} \mye 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 {\sl{ottimale}} rispetto al tempo, al numero di operazioni, alla stabilit\mya e all'accuratezza. Ovviamente ad un maggior numero di operazioni non \mye detto che corrisponda una miglior accuratezza dei risultati. La complessit\mya computazionale conta il numero di operazioni necessarie per fornire l'approssimazione desiderata.

\subsection{Esempio: algoritmo di Horner.}
Ci poniamo il problema di valutare il polinomio
\begin{equation}\label{f1}
p(x)=a_0+a_1 \cdot x+ \ldots + a_n \cdot x^n
\end{equation}
in un punto $x$. 

Osserviamo che 
\begin{equation}\label{f2}
p(x)= a_0+x\cdot(a_1+x\cdot(a_2+\ldots+x\cdot(a_{n-1}+x\cdot a_n)))
\end{equation}

Supponiamo sia $a=(a_0,\ldots,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 ({\ref{f1}}), il secondo che effettua la stessa operazione come descritto in ({\ref{f2}}) calcolando dapprima $s_1=a_{n-1} + x \cdot a_n$, poi  $s_2=a_{n-2} + x \cdot s_1$ e cos\myi via. In Matlab avremo allora 

\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

function s=algoritmo1(a,x)

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

\end{Verbatim}
e
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}

Se lanciamo il codice {\tt{demo\_horner}} per la valutazione di $p(x)=1+2 \cdot x + 3 \cdot x^2 +4 \cdot x^3$ in $x=\pi$
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

clear all;

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


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

format long;
y1
y2

\end{Verbatim}
otteniamo
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}

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


\section{Facoltativo: Un esempio sulla soluzione delle equazioni di secondo grado}

Consideriamo l'equazione di secondo grado
\begin{equation}\label{eq.st.1}
x^2 - 2\sqrt{3} x + 3 = 0
\end{equation}

\noindent
Si nota subito che 
\begin{equation}\label{eq.st.1}
x^2 - 2\sqrt{3} x + 3 = (x-\sqrt{3})^2
\end{equation}

\noindent
e quindi (\ref{eq.st.1}) ha un'unica soluzione $\sqrt{3}$ con molteplicit\mya 2, poich\mye la derivata di $x^2 - 2\sqrt{3} x + 3$ \mye la funzione  $2 x - 2\sqrt{3} x$ si annulla in $\sqrt{3}$ (cf. {\cite[p.420]{COM}}).

Per quanto riguarda il condizionamento di una radice $\alpha$ (cf. {\cite[p.420]{COM}}), si prova che se $\alpha({\epsilon})$ \mye la radice di 
\begin{equation}\label{eq.st.2}
P_{\epsilon}(z)=P(z)+\epsilon g(z), \,\,\, \epsilon>0
\end{equation}
allora se $\alpha$ \mye semplice
\begin{equation}\label{eq.st.3}
\alpha({\epsilon}) \approx \alpha + \epsilon \left(    \frac{ -g(\alpha)}{P^1(\alpha)}\right)
\end{equation}
mentre se $\alpha$ \mye multipla con moltiplicit\mya $m$, cio\mye 
$$P(\alpha)= \ldots=P^{(m-1)}(\alpha)=0 \,\,\, P^{(m)}(\alpha) \neq 0$$
allora esiste una radice $\alpha({\epsilon})$ di $P_{\epsilon}$ tale che
\begin{equation}\label{eq.st.4}
\alpha({\epsilon}) \approx \alpha + \epsilon^{1/m} \left(    \frac{ -m!g(\alpha)}{P^{(m)}(\alpha)}\right)
\end{equation}

Approssimiamo in (\ref{eq.st.1}) il valore $\sqrt{3}$ con $1.7321$ fornito da Matlab usando il cosidetto {\tt{format short}}. In questo caso, posto $g(z)=2z$, $$\epsilon=-(\sqrt{3}-1.7321) \approx 4.9192 \cdot 10^{-5} \approx 5 \cdot 10^{-5}$$ abbiamo che 
\begin{equation}\label{eq.st.5}
\alpha({\epsilon}) \approx \alpha + \epsilon^{1/m} \left(    \frac{ -m!\alpha}{P^{(2)}(\alpha)}\right)
\end{equation}
\noindent
ed essendo $P^{(2)}(z)=2$, $m=2$ (molteplicit\mya della radice $\alpha=\sqrt{3}$), ricaviamo
\begin{equation}\label{eq.st.5}
\alpha({\epsilon}) \approx \alpha + (5 \cdot 10^{-5})^{1/2} \left(    \frac{ -2! 2 \sqrt{3}}{2}\right)
\end{equation}
in cui 
$$|(5 \cdot 10^{-5})^{1/2} \left(    \frac{ -2! 2\sqrt{3}}{2}\right)^{1/2} |\approx 0.02449489742783$$
Quindi ci si aspetta che una delle radici di 
$$P_{\epsilon}(z)=x^2 - 2 \cdot 1.7321 \cdot x + 3$$
disti circa $0.02449489742783$ da $\sqrt{3}$.

In effetti, utilizzando semplici modifiche del programma {\tt{radicisecondogrado}}, si hanno due radici,
$$x_1 \approx 1.7451541181241765000, \,\,\, x_2 \approx 1.7190458818758234000.$$
ed \mye, come si vede facilmente dalla shell di Matlab/Octave:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}] 

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

\end{Verbatim}


\section{Esercizio}
Date le funzioni 
$$f_1(x)=1-\sqrt{1-x^2}$$
$$f_2(x)=1-x$$
si calcoli analiticamente il {\sl{condizionamento}}
$${\cal{K}}_1(x)=\frac{|  x \cdot f_1^{\prime}(x)|}{|f_1(x)|} $$
$${\cal{K}}_2(x)=\frac{|  x \cdot f_2^{\prime}(x)|}{|f_2(x)|} $$
Utilizzando Matlab/Octave, si plottino in scala semilogaritmica  i grafici di ${\cal{K}}_1$ e ${\cal{K}}_2$ nel set di punti $x=-1+10^{(-3)},-1+2\cdot10^{(-3)},\ldots,1-2 \cdot 10^{(-3)},1-10^{(-3)},$. Si pu\myo dire da questi dove la funzione \mye malcondizionata (cio\mye assume valori di condizionamento {\sl{alti}})?

\section{Alcuni esercizi}
{\smsp}
\begin{enumerate}

\item {\bf{Esercizio facoltativo di media difficolt\mya}}. Esistono pi\myu successioni che approssimano $\pi$ e sono attribuite ad Archimede. Ad esempio, osservato che questi non \mye altro che l'area del cerchio avente raggio $1$, si inscrivono nel cerchio al variare di $p=2,3,\ldots$ dei poligoni regolari aventi $2^p$ lati. Osserviamo che per $p=2$, si deve determinare l'area del quadrato inscritto, che con facili conti \mye 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).

 \begin{figure}\label{fig1s}
 \centering
   {\includegraphics[scale=0.35,clip]{archimede_8.eps}}
   {\includegraphics[scale=0.35,clip]{archimede_16.eps}}
 \caption{Grafico che illustra le suddivisioni considerate dall'algoritmo di Archimede, per $p=3$, $p=4$.}
 \end{figure}

Se indichiamo con $\theta$ l'angolo al centro di ogni triangolo, si verifica che la sua area \mye $\frac{\sin {(\theta)}}{2}$. In particolare per un poligono di $2^p$ lati, si ha $\theta_p=\frac{2\pi}{2^p}$ e l'area del poligono inscritto, costituito da $2^p$ triangoli uguali, \mye 
$$I_p=2^p \, \frac{\sin {(\theta_p)}}{2} = 2^p \, \frac{\sin {(\frac{2\pi}{2^p})}}{2}.  $$
In realt\mya Archimede us\myo questa ed altre idee per cercare di approssimare $\pi$ con stime dall'alto e dal basso, via anche poligoni circoscritti. Come citato in {\cite{WIKAR}}:

{\sl{Nel breve lavoro {\rm{La misura del cerchio}} viene dimostrato anzitutto che un cerchio \mye equivalente a un triangolo con base eguale alla circonferenza e altezza eguale al raggio. Tale risultato \mye 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\myo approssimare arbitrariamente il rapporto tra circonferenza e diametro di un cerchio dato, rapporto che oggi si indica con $\pi$. Le stime esplicitamente ottenute limitano questo valore fra 22/7 e 3 + 10/71}}.

Secondo {\cite{WOL}}, Archimede us\myo un altro algoritmo per approssimare $\pi$. Visto che $2\pi $ non \mye altro che la lunghezza della circonferenza del cerchio avente raggio uguale a $1$ e questa \mye maggiore del perimetro $a_k$ di un qualsiasi poligono regolare di $n=6 \cdot 2^k$ lati inscritto e minore del perimetro di un qualsiasi poligono regolare di $n=6 \cdot 2^k$ lati circoscritto $b_k$, misurando $a_k$ e $b_k$ si pu\myo affermare che $$a_k \leq 2\pi \leq b_k. $$

\noindent
Si pu\myo inoltre provare che 
$$a_k=2n \,  \sin {\left(\frac{\pi}{n} \right)},$$
$$b_k=2n \, \tan {\left(\frac{\pi}{n}\right)}.$$
Stimare $\pi$ come $$a_k \leq 2 \, \pi \leq b_k.$$ 

\noindent
Per quale $n$ si riesce a determinare $\pi$ con $10$ cifre decimali esatte? Pensare se sia giusto dire $b_k-a_k < 10^{-10}$ allora $\pi$ \mye approssimato da $a_k$ con 10 cifre decimali esatte.

 \begin{figure}\label{fig_archimede}
 \centering
   {\includegraphics[scale=10.0,clip]{stamp_italian_small_gimp.eps}}
 \caption{Un francobollo raffigurante Archimede (287ac-212ac).}
 \end{figure}

\noindent
{\bf{A.}} Si approssimi $\pi$ 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.

\noindent
{\bf{B.}} Stimare $\pi$ come $$a_k \leq \pi \leq 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 $\pi$ con $10$ cifre decimali esatte? Per $k=0,\ldots,4$ si ha
\begin{eqnarray}
3.00000 & \leq \pi \leq & 3.46410 \\	
3.10583 & \leq \pi \leq & 3.21539 \\
3.13263 & \leq \pi \leq & 3.15966 	\\
3.13935 & \leq \pi \leq & 3.14609 	\\
3.14103 & \leq \pi \leq & 3.14271
\end{eqnarray}

\noindent
{\bf{Osservazione}}: il fatto che ci sia un $\pi$ nei secondi membri di alcune espressioni non \mye {\sl{un gatto che si morde la coda}}. Infatti \mye solo una trascrizione in radianti di un angolo che ha un significato geometrico indipendente dalla quantit\mya $\pi$. In questi esercizi, per semplificare la questione si usa per\myo $\pi$ per calcolare $\pi$, cosa che dimostra una volta in pi\myu l'abilit\mya di Archimede, che calcol\myo tali quantit\mya con metodi elementari e geometrici.

Si cita il fatto che per $k=4$ Archimede riusci' ad affermare che 
$$\frac{223}{71} \leq \pi \leq \frac{22}{7}.$$
{\smsp}
% ESERCIZIO 5.
\item {\bf{Esercizio, facile facoltativo}}. Posto $h=10^{-1},10^{-2},\ldots,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$?
{\smsp}
% ESERCIZIO 6.
\item {\bf{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.
{\smsp}
% ESERCIZIO 7.
\item {\bf{Esercizio, facile facoltativo}}. Fissato $\theta=0:(2\pi / 1001):2 \pi$, si valuti $\sin^2{(\theta)} + \cos^2{(\theta)}$ e si calcoli l'errore assoluto rispetto il valore $1$.
\end{enumerate}
























































D'altro canto, spesso la funzione $f$ non \mye valutabile esattamente ma approssimata da una certa $g$. Il valore 
$$\frac{g({x})-f(x)}{f(x)}$$
si chiama  {\sl{errore relativo analitico o di troncamento}}.

Una volta definita $g$, il suo calcolo viene effettuato da un certo algoritmo ottenendo una approssimazione ${\overline{g}}$. Per $g(x) \neq 0$, la quantit\mya 
 $$\frac{{\overline{g}}({x})-{\overline{g}}(x)}{g(x)}$$
si chiama  {\sl{errore relativo algoritmico}}.

Da questa analisi, si capisce che usualmente, dovendo calcolare $f(x)$ si ottiene, per varie ragioni, la quantit\mya ${\overline{g}}({\overline{x}})$.

\section{Condizionamento e stabilit\mya}

Supponiamo che la funzione scalare $f$ sia valutabile esattamente.  La sua valutazione in un punto $x$ si dice {\sl{bencondizionata}} se a piccole variazioni relative sui dati corrispondono piccole variazioni sui risultati. In altre parole, se $x \neq 0$ e $f(x) \neq 0$ si ha 
  $$\frac{|f({\overline{x}})-f(x)|}{|f(x)|} \approx \frac{|({\overline{x}})-(x)|}{|(x)|}.$$
Osserviamo che abbiamo supposto che $f$ sia valutabile esattamente e quindi il bencondizionamento dipende esclusivamente da $f$ e dal dato $x$ e non dall'algoritmo utilizzato.
Sottolineiamo che la valutazione di $f$ pu\myo essere malcondizionata per certi dati $x$ e bencondizionata per altri $x$.

Un algoritmo per la valutazione di $f$ si dice {\sl{stabile}} se amplifica {\sl{poco}} gli errori di arrotondamento introdotti dalle singole operazioni. Osserviamo che la valutazione di una funzione $f$ bencondizionata, pu\myo dare risultati non sufficientemente corretti se effettuata da un algoritmo instabile. 


\section{Facoltativo. Alcuni esercizi}
{\smsp}
\begin{enumerate}

% ESERCIZIO 1

\item {\bf{Esercizio facile}}. Si calcoli l'errore relativo tra $1$ e il valore che si ottiene
valutando in Matlab/Octave
$${\frac {(1+\eta)-1} {\eta} }$$   
con $\eta=10^{-1}, 10^{-2}, \ldots, 10^{-15}$. Si consideri poi 
$\eta=8.8817841970012523E-16$ e si calcoli la medesima quantit\'a, 
giustificando i risultati ottenuti {\cite[p. 5, problema 3]{FRO}}. 

\noindent
Suggerimento: quanto vale {\tt{eps}}?
{\smsp}
% ESERCIZIO 2

\item {\bf{Esercizio facile}}. Siano $f:=\tan$ e $g:=\arctan$. Si consideri la funzione composta $$f 
\cdot 
g(x):=f(g(x))=\tan(\arctan(x))=x, \,\,x \in (-\infty,+\infty).$$ Si 
calcolino
$$x=10^{k}, {\mbox{ per }}  k=-20+40h, {\mbox{ e }} 
h=0,0.01,0.02,\ldots,1$$ 
e si valuti l'errore relativo compiuto. Si pu\'o dire che anche 
{\sl{numericamente}} $\tan(\arctan(x))=x$?
{\smsp}
% ESERCIZIO 4

\item {\bf{Esercizio facoltativo di media difficolt\mya}}. Esistono pi\myu successioni che approssimano $\pi$ e sono attribuite ad Archimede. Ad esempio, osservato che questi non \mye altro che l'area del cerchio avente raggio $1$, si inscrivono nel cerchio al variare di $p=2,3,\ldots$ dei poligoni regolari aventi $2^p$ lati. Osserviamo che per $p=2$, si deve determinare l'area del quadrato inscritto, che con facili conti \mye 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).

 \begin{figure}\label{fig1s}
 \centering
   {\includegraphics[scale=0.35,clip]{archimede_8.eps}}
   {\includegraphics[scale=0.35,clip]{archimede_16.eps}}
 \caption{Grafico che illustra le suddivisioni considerate dall'algoritmo di Archimede, per $p=3$, $p=4$.}
 \end{figure}

Se indichiamo con $\theta$ l'angolo al centro di ogni triangolo, si verifica che la sua area \mye $\frac{\sin {(\theta)}}{2}$. In particolare per un poligono di $2^p$ lati, si ha $\theta_p=\frac{2\pi}{2^p}$ e l'area del poligono inscritto, costituito da $2^p$ triangoli uguali, \mye 
$$I_p=2^p \, \frac{\sin {(\theta_p)}}{2} = 2^p \, \frac{\sin {(\frac{2\pi}{2^p})}}{2}.  $$
In realt\mya Archimede us\myo questa ed altre idee per cercare di approssimare $\pi$ con stime dall'alto e dal basso, via anche poligoni circoscritti. Come citato in {\cite{WIKAR}}:

{\sl{Nel breve lavoro {\rm{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 \mye 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\myo approssimare arbitrariamente il rapporto tra circonferenza e diametro di un cerchio dato, rapporto che oggi si indica con $\pi$. Le stime esplicitamente ottenute limitano questo valore fra 22/7 e 3 + 10/71}}.

Secondo {\cite{WOL}}, Archimede us\myo un altro algoritmo per approssimare $\pi$. Visto che $2\pi $ non \mye altro che la lunghezza della circonferenza del cerchio avente raggio uguale a $1$ e questa \mye maggiore del perimetro $a_k$ di un qualsiasi poligono regolare di $n=6 \cdot 2^k$ lati inscritto e minore del perimetro di un qualsiasi poligono regolare di $n=6 \cdot 2^k$ lati circoscritto $b_k$, misurando $a_k$ e $b_k$ si pu\myo affermare che $$a_k \leq 2\pi \leq b_k. $$

\noindent
Si pu\myo inoltre provare che 
$$a_k=2n \,  \sin {\left(\frac{\pi}{n} \right)},$$
$$b_k=2n \, \tan {\left(\frac{\pi}{n}\right)}.$$
Stimare $\pi$ come $$a_k \leq 2 \, \pi \leq b_k.$$ 

\noindent
Per quale $n$ si riesce a determinare $\pi$ con $10$ cifre decimali esatte? Pensare se sia giusto dire $b_k-a_k < 10^{-10}$ allora $\pi$ \mye approssimato da $a_k$ con 10 cifre decimali esatte.

 \begin{figure}\label{fig_archimede}
 \centering
   {\includegraphics[scale=10.0,clip]{stamp_italian_small_gimp.eps}}
 \caption{Un francobollo raffigurante Archimede (287ac-212ac).}
 \end{figure}

\noindent
{\bf{A.}} Si approssimi $\pi$ 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.

\noindent
{\bf{B.}} Stimare $\pi$ come $$a_k \leq \pi \leq 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 $\pi$ con $10$ cifre decimali esatte? Per $k=0,\ldots,4$ si ha
\begin{eqnarray}
3.00000 & \leq \pi \leq & 3.46410 \\	
3.10583 & \leq \pi \leq & 3.21539 \\
3.13263 & \leq \pi \leq & 3.15966 	\\
3.13935 & \leq \pi \leq & 3.14609 	\\
3.14103 & \leq \pi \leq & 3.14271
\end{eqnarray}

\noindent
{\bf{Osservazione}}: il fatto che ci sia un $\pi$ nei secondi membri di alcune espressioni non \mye {\sl{un gatto che si morde la coda}}. Infatti \mye solo una trascrizione in radianti di un angolo che ha un significato geometrico indipendente dalla quantit\mya $\pi$. In questi esercizi, per semplificare la questione si usa per\myo $\pi$ per calcolare $\pi$, cosa che dimostra una volta in pi\myu l'abilit\mya di Archimede, che calcol\myo tali quantit\mya con metodi elementari e geometrici.

Si cita il fatto che per $k=4$ Archimede riusci' ad affermare che 
$$\frac{223}{71} \leq \pi \leq \frac{22}{7}.$$
{\smsp}
% ESERCIZIO 5.
\item {\bf{Esercizio, facile facoltativo}}. Posto $h=10^{-1},10^{-2},\ldots,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$?
{\smsp}
% ESERCIZIO 6.
\item {\bf{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.
{\smsp}
% ESERCIZIO 7.
\item {\bf{Esercizio, facile facoltativo}}. Fissato $\theta=0:(2\pi / 1001):2 \pi$, si valuti $\sin^2{(\theta)} + \cos^2{(\theta)}$ e si calcoli l'errore assoluto rispetto il valore $1$.
\end{enumerate}






\begin{center}
\begin{table}\label{tabden}
\small
\begin{tabular}{| c | c | c | c |}
\hline 
{\mbox{precisione}} & {\mbox{denorm.}} & {\mbox{norm.}} & {\mbox{in decimali}} \\
 {\mbox{singola}} & $\pm 2^{-149} \rightarrow (1-2^{-23}) \cdot 2^{+126}$ &  $\pm 2^{-126} \rightarrow (1-2^{-23}) \cdot 2^{127}$  & $$  $\pm \approx 10^{-44.85} \rightarrow \approx 10^{38.53}$  \\
{\mbox{doppia}} &  $\pm 2^{-1074} \rightarrow (1-2^{-52}) \cdot 2^{+1022}$ &  $\pm 2^{-1022} \rightarrow(2-2^{-52}) \cdot 2^{1023}$  & $\pm \approx 10^{-323.3} \rightarrow \approx 10^{308.3}$  \\
\hline
\end{tabular}
\caption{Tabella rappresentante il range dei numeri denormalizzati e normalizzati}
\end{table}
\end{center}

\end{remark}








\subsection{Facoltativo (argomento difficile). Consistenza, stabilit\mya e convergenza}
Vediamo i dettagli del concetto di stabilit\mya e li legheremo ad altre due propriet\myas, quelle di {\sl{consistenza}} e {\sl{convergenza}} (cf.{\cite[p.39-41]{MATNUM}}).
Supponiamo di voler risolvere il problema 
\begin{equation}\label{stprob}
F(x,d)=0, 
\end{equation} 
e che questo sia ben posto, cio\mye a piccole variazioni di $d$ corrispondano piccole variazioni di $x$. Nella precedente scrittura, la variabile $d$ \mye fissata e si cerca $x$ che risolva il problema (\ref{stprob}).

Un metodo numerico per la soluzione approssimata di $F(x,d)$ costruisce una successione di problemi $$F_n(x_n,d_n)=0, \hsp n \geq 1$$
dipendenti da $n$, col desiderio che sia $$x_n \rightarrow x, \hsp n \rightarrow \infty.$$
Se il dato $d$ \mye ammissibile per $F_n$, cio\mye ha senso risolvere $F_n(x,d)=0$, si dice che il metodo \mye {\sl{consistente}} se per $x$ soluzione del problema $F(x,d)=0$ si ha
$$F_n(x,d)=F_n(x,d)-F(x,d) \rightarrow 0, \hsp n \rightarrow \infty.$$
Il metodo si dice {\sl{fortemente consistente}} se in particolare 
$$F_n(x,d)=F_n(x,d)-F(x,d) = 0, \hsp \forall n.$$

Come detto in precedenza, un metodo \mye stabile (talvolta si dice {\sl{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\mya dai dati cio\mye
$$\forall \eta > 0, \hsp \exists K_n(\eta,d_n): \hsp \|\delta d_n\| < \eta \Rightarrow  \|\delta x_n\| \leq K_n(\eta,d_n)\|\delta d_n\|.$$

Il proposito di un metodo 	\mye ovviamente di costruire, attraverso problemi numerici$$F_n(x_n,d_n)=0,$$ una successione di $x_n$ che si avvicini sempre pi\myu a $x$ al crescere di $n$, cio\mye il metodo sia {\sl{convergente}}. In termini pi\myu rigorosi, se $F(x(d),d)=0$ e $F_n(x_n(d+\delta d_n),d+\delta d_n)=0$
\begin{eqnarray}
\forall \epsilon &> 0& \hsp  \exists n_0(\epsilon), \exists \delta(n_0,\epsilon) > 0: \forall n > n_0(\epsilon), \forall \|\delta d_n\| < \delta(n_0,\epsilon) \nonumber \\ &\hsp& \Rightarrow \|x(d) -x_n(d+\delta d_n)\| \leq \epsilon.
\end{eqnarray}

\noindent Il teorema di Lax-Richtmyer mostra come queste tre propriet\mya siano interconnesse, cio\mye 
\vspace{0.25cm}
\begin{theorem}
Per un metodo numerico consistente, la stabilit\mya \mye equivalente alla convergenza.
\end{theorem}
\vspace{0.25cm}

