
% Copyright (C) 2012 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{geometry}
\usepackage[colorlinks=true,linkcolor=purple,urlcolor=purple]{hyperref}
\usepackage{float}
\floatstyle{ruled}
\newfloat{program}{thp}{lop}
\floatname{program}{Codice}
\pdfpagewidth 8.5in
\pdfpageheight 11in
\geometry{left=3.5cm, right = 3.5cm, top = 3cm, bottom = 3cm}
% \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{\myed}{\'o }
\newcommand{\myu}{\`u }
\newcommand{\myur}{\'u }
\newcommand{\mya}{\`a }
\newcommand{\myo}{\`o }
\newcommand{\myod}{\'o }
\newcommand{\myi}{\`i }
\newcommand{\myid}{\'i }
\newcommand{\myfl}{\mbox{fl}}
\newcommand{\smsp}{\vspace{0.2cm}}
\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
  }}
\renewcommand{\sectionmark}[1]{\markright{\thesection\#1}}
\fancyhf{}
\fancyhead[LE]{\bf\nouppercase{\thepage \hfill  Laurea in Matematica A. A. 2012--2013}}
\fancyhead[RO]{\bf\nouppercase{Laboratorio di Calcolo Numerico \hfill   \thepage }}
\pagestyle{fancy}


\title{\bf Sulla propagazione degli errori \thanks{Ultima revisione: \today}}
\author{}
\author{{A. Sommariva
{\thanks{Dipartimento di Matematica, Universit\`a degli Studi di  Padova, via Trieste 63, 35121 Padova, Italia ({\tt alvise@math.unipd.it})}}}  {    A. Mart\'{\i}nez
{\thanks{Dipartimento di Matematica, Universit\`a degli Studi di Padova,
 via Trieste 63, 35121 Padova, Italia ({\tt angeles.martinez@unipd.com})}} }
}
\begin{document}
\maketitle
\section{Esercitazione sulla propagazione degli errori}
\subsection{Facoltativo: Rappresentazione dei numeri al calcolatore}

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

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

Si intende per {\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
occupano una {\sl{parola}} da 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{1}
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}}. La scrittura $1.F$ significa $1+F$ dove $F=\sum^{+\infty}_{k=1}a_k 2^{-k}$, con $a_k=0$ oppure $a_k=1$. Si noti che delle possibili infinite cifre decimali del numero, solo un numero limitato di esse (Fbits cifre) vengono memmorizzate all'interno del calcolatore. In particolare per la singola e la doppia precisione si hanno le seguenti quantit\mya
%Alcuni esempi
{\smsp}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
Precisione & Ebits & Bias & Fbits \\
Singola    & 8         & 127  & 23        \\
Doppia     & 11        & 1023 & 52         \\
\hline
\end{tabular}
\end{center}
{\smsp}
\noindent
Per capire meglio tale notazione consideriamo il link

\vspace{0.2cm}
\begin{center}
\href{http://babbage.cs.qc.edu/IEEE-754/Decimal.html}{http://babbage.cs.qc.edu/IEEE-754/Decimal.html}
\end {center}
\vspace{0.2cm}

\noindent
in cui viene convertito un numero decimale in forma binaria, secondo lo standard IEEE 754.
Partiamo con un esempio semplice.

\noindent
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*2^0+1*2^1+0*2^2+0*2^3+0*2^4+0*2^5+0*2^6+1*2^7=2+128=130$,
l'esponente vale $130-127=3$. Quindi $2^{E-Bias}=2^3=8$.
{\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*2^{-1}+0*2^{-2}+\ldots=0.5$. Ora
($1.F$) si legge come $1+F=1.5$.
\end{enumerate}
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*2^0+0*2^1+1*2^2+0*2^3+0*2^4+0*2^5+0*2^6+1*2^7=1+4+128=133$,
l'esponente vale $133-127=6$. Quindi $2^{E-Bias}=2^6=64$. Viene da chiedersi perch\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*2^{-1}+1*2^{-2}+
1*2^{-3}+1*2^{-4}+1*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}
In altre parole la rappresentazione di $126$ \mye esatta ed il numero \mye
rappresentato come $126=64\,\cdot\,1.96875$.

Dalla tabella riguardante la precisione singola e i bit, si deduce che i
numeri in virgola mobile non sono dunque un insieme denso, come sono
invece i numeri reali, bens\myi un insieme discreto di valori sull'asse
reale, in posizioni non equispaziate.  Infatti, la distanza tra un numero
in virgola mobile $x$ ed il successivo numero macchina $x_+$, \mye:
$$|x - x_+| = 2^{\mbox{-nbitsF}}  \,\cdot\,  2^{\mbox{E-bias}}.$$

Si dice precisione di macchina, la met\mya della distanza
tra $1$ ed il successivo numero macchina.
In Matlab tale distanza \mye rappresentata da ${\tt{eps}}$ ed \mye pari a $2^{-52}$. 

La cardinalit\mya dell'insieme di numeri macchina \mye piuttosto elevata: in dop\-pia pre\-ci\-sio\-ne i nu\-meri
macchina sono $$2^{64} \approx 1.844674407370955e+019.$$
%Vediamo di
%calcolarlo. Si scriva sulla shell di Linux
%\begin{center}
%\tt{pico myeps.m}
%\end{center}
%e nel file il codice Matlab
%\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]
%
%myeps=1;
%while 1+myeps > 1;
%      myeps=myeps/2;
%end
%myeps=2*myeps;
%%format long e
%myeps-eps

%\end{Verbatim}
%Il codice Matlab parte con {\tt{myeps=1}} e continua a dividere per $2$
%finch\'e trova un numero $x$ tale che sommato a $1$ d\mya (numericamente,
%non analiticamente) $1$. Vista la definizione di {\tt{eps}} alla fine bisogna porre
%${\mbox{myeps}}=2*{\mbox{myeps}}$ (altrimenti ${\mbox{myeps}}<{\mbox{eps}}$). E infatti ${\mbox{myeps-eps}}=0$.
%
\subsection{Facoltativo: Su {\tt{eps}}, {\tt{realmax}}, {\tt{realmin}}, overflow e underflow}
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-{\mbox{bias}} } $ ed in Matlab \mye
espresso dalla variabile {\tt{realmin}}.
Il limite superiore \mye pari a $2^{2^{\mbox{Ebits}}-2-bias}*(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$,
{\sl{Ebits}}$=11$, ed {\sl{Fbits}}$=52$. Conseguentemente
\begin{eqnarray}
{\mbox{realmax}} & = & 2^{2^{\mbox{Ebits}}-2-bias} \, \cdot \, (1+(1- 2^{\mbox{-Fbits}} )) \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*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*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 .* (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}
Lo standard richiede che il risultato delle operazioni di addizione,
sottrazione, moltiplicazione e divisione sia arrotondato esattamente,
cio\mye il risultato deve essere calcolato esattamente e poi arrotondato
al numero in virgola mobile pi\myu vicino.

Problemi in tal senso provengono dalla sottrazione. Per questo motivo, nei
microprocessori moderni la sottrazione viene effettuata utilizzando la
tecnica bit di guardia. Non tratteremo questo punto dettagliatamente. Per
curiosit\mya si consulti il sito
{\smsp}
\begin{center}
\href{http://babbage.cs.qc.edu/IEEE-754/References.xhtml}{http://babbage.cs.qc.edu/IEEE-754/References.xhtml}
\end{center}
{\smsp}
{\subsection{Facoltativo: Propagazione degli errori di arrotondamento}}

Ogni calcolatore pu\myo rappresentare esattamente solo un numero finito di numeri reali, i cosidetti numeri
macchina. Per i rimanenti, seguendo uno standard quale IEEE 754, pu\myo fornirne solo un'ap\-pros\-si\-ma\-zio\-ne. Di conseguenza, le operazioni aritmetiche di
base {\sl{possono}} in generale essere soggette ad un errore nel risultato ed \mye dunque necessario conoscere l'entit\mya di quest'ultimo e l'effetto
che pu\myo avere in un algoritmo.

\noindent
Questo problema non \mye da trascurarsi. Si possono trovare in internet vari siti in cui queste approssimazioni, dovute ad esempio al cambiamento di
unit\mya di misura, hanno portato a disastri quali la perdita di satelliti o la distruzione di stazioni petrolifere come indicato ad esempio in
{\cite{HUC}}.

\noindent
Passiamo quindi in dettaglio all'analisi del problema. Sia $x \neq 0$ un numero reale, che in virgola mobile normalizzata si pu\myo scrivere come 
$$x = (-1)^S \cdot 2^{(E-{\mbox{bias}})} \cdot (1.F + \delta)$$
in cui supponiamo  $|\delta| < 2^{(-{\mbox{Fbits}})} / 2$. 
%Ricordiamo che  la scrittura $1.F$ significa $1+F$ dove $F=\sum^{+\infty}_{k=1}a_k 2^{-k}$, con $a_k=0$ oppure $a_k=1$.
Se si considera precisione singola, secondo lo standard IEEE 754 si avr\mya ${\mbox{Ebits}}=8$, ${\mbox{bias}}=127$, ${\mbox{Fbits}}=23$. 
Allora, la rappresentazione di $x$ all'interno del calcolatore, sar\mya
$${\mbox{fl}}(x)  =   (-1)^S  \cdot  2^{(E-{\mbox{bias}})} \cdot  1.F$$
con $F=\sum^{Fbits}_{k=1}a_k 2^{-k}$ e dunque con un errore relativo (di arrotondamento)

$$ {\frac{|x-{\mbox{fl}}(x)|}{|x|}} = {\frac{| 2^{(E-{\mbox{bias}})} \cdot \delta |}{|2^{(E-{\mbox{bias}})} \cdot (1.F + \delta)|}}
= {\frac{| \delta |}{|1.F + \delta|}} $$

\noindent
Ricordiamo che 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.


\noindent
Per quanto riguarda le operazioni aritmetiche fondamentali si pu\myo dimostrare che la somma, la divisione e la moltiplicazione producono un errore
relativo piccolo, mentre la
sottrazione pu\myo anche produrre un errore relativo grande rispetto al risultato (ci\myo avviene quando i due operandi sono molto vicini tra di loro
e si ha dunque una perdita di cifre significative nel risultato).

\noindent
Pi\myu precisamente se ${\mbox{fl}}(x)$ \mye il numero macchina che {\sl{corrisponde}} a $x$, denotati 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$, $:$, posto
\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}

\noindent
si ha dopo qualche non banale conto (cf. {\cite[p.78]{COM}})
\begin{eqnarray}
\epsilon_{x,y}^{\oplus} &\approx&  \mylbar {\frac{x}{x+y} } \myrbar \epsilon_x + \mylbar {\frac{y}{x+y} } \myrbar \epsilon_y \\
\epsilon_{x,y}^{\ominus} &\approx&  \mylbar {\frac{x}{x-y} } \myrbar \epsilon_x + \mylbar {\frac{y}{x-y} } \myrbar \epsilon_y \\
\epsilon_{x,y}^{\otimes} &\approx& \epsilon_x + \epsilon_y \\
\epsilon_{x,y}^{\oslash} &\approx& \mylbar \epsilon_x - \epsilon_y \myrbar
\end{eqnarray}

\noindent
il che mostra il pericolo di perdita di accuratezza nella somma e nella sottrazione qualora rispettivamente $x+y \approx 0$ e $x-y \approx 0$.



\noindent

\section{Calcolo di una radice di una equazione di secondo grado}

Vediamo un esempio concreto in cui l'introdurre una sottrazione potenzialmente pericolosa conduce effettivamente a problemi di
innacuratezza 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$.
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 stabile}} (o {\sl{algoritmo instabile}}) {\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}

\noindent
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}$.
Secondo {\cite{COMESE}}, p. 11 \mye notevole l'esperimento in cui $$p=4.999999999995 \, 10^{+4}, \, q=10^{-2}$$ avente soluzione esatta $10^{-7}$.
Si osservi che in entrambi i casi effettivamente $p \gg q$.

\noindent
Passiamo ora all'implementazione e verifica numerica di questi due tests, osservando che i problemi relativi sono comunque ben condizionati.

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

% p=4.999999999995*10^(+4); q=10^(-2); sol=10^(-7);
p=1000; q=0.018000000081; sol=0.9*10^(-5);

% ALGORITMO 1
s=p^2;
t=s+q;
if t >=0
   u=sqrt(t);
else
   fprintf('\n \t [RADICI COMPLESSE]');
end
s1=-p+u;

% ALGORITMO 2
s=p^2;
t=s+q;
if t >=0
   u=sqrt(t);
else
   fprintf('\n \t [RADICI COMPLESSE]');
end
v=p+u;
t1=q/v;

fprintf('\n \t [ALG.1] [1]: %10.19f',s1);
fprintf('\n \t [ALG.2] [1]: %10.19f',t1);
if length(sol) > 0 & (sol ~= 0)
    relerrs1=abs(s1-sol)/abs(sol);
    relerrt1=abs(t1-sol)/abs(sol);
    fprintf('\n \t [REL.ERR.][ALG.1]: %2.2e',relerrs1);
    fprintf('\n \t [REL.ERR.][ALG.2]: %2.2e',relerrt1);
end

\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.3cm}}

\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*10^(+4); q=10^(-2); sol=10^(-7);

\end{Verbatim}

{\vspace{0.3cm}}

\noindent
{\bf{Esercizio}}. 
Si determini una formula stabile alternativa per il calcolo delle radici dell'equazione di secondo grado $x^2 + 2 p x -q$. \\

\noindent {Suggerimento}: si consideri la relazione che esiste tra il prodotto delle radici e il coefficiente $q$.

\section{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}}.

\subsection{Successioni convergenti a $\pi$: tests numerici}

\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$).
Non \mye difficile osservare che 
$$\left\{ \begin{array} {l}
y_1=1, \, \, y_2=2 \\
y_{n+1} = 
{ 
\frac 
{2^{{ \frac {1} {2} } }  |y_n|} 
{ \sqrt{1 + \sqrt{1-4^{1-n} \cdot y_n^2}}
 } 
}
\end{array} \right.$$

\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.2cm}}

\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.3cm}}

{\subsubsection{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 con 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}

\section{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} \approx 0.3679$ e pi\myu in generale integrando per parti, cio\mye osservando che
$$\int_a^b f^{\prime}(x) \, g(x) dx = f(x) \, g(x)\mid_a^b - \int_a^b f(x) \, g^{\prime}(x)  dx, $$
ricaviamo per $f(x)=\exp{(x)}$ e $g(x)=x^{n+1}$ e dalla definizione di $I_n$
\begin{eqnarray}\label{322}
I_{n+1}&=&\exp{(-1)} \int_0^1 x^{n+1} \, \exp{(x)} \, dx \nonumber \\
&=& \exp{(-1)} \left( x^{n+1} \, \exp{(x)}  \mid_0^1 - (n+1) \int_0^1 x^n \, \exp{(x)} 
\, dx  \right) \nonumber\\ 
&=& \exp{(-1)}  \left( 1^{n+1} \cdot \exp{(+1)} - 0 \cdot \exp{(0)}  - (n+1) \, \exp{(+1)}  \cdot \exp{(-1)}  \int_0^1 x^n \, \exp{(x)} \, dx  \right) \nonumber\\ 
&=& 1- (n+1) { \exp{(-1)}  \int_0^1 x^n \, \exp{(x)} \, dx}  \nonumber\\ 
&=&1-(n+1)\, I_n.
\end{eqnarray}

\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.3cm}}

{\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.3cm}}
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

% SUCCESSIONE RICORRENTE.

% SUCCESSIONE "s_n" (IN AVANTI).
s(1)=exp(-1);   
for n=1:99
    s(n+1)=1-(n+1)*s(n);
end

% SUCCESSIONE "t_n" (ALL'INDIETRO).
m=500; M=2*m; % PER CALCOLARE I VALORI TRA 
              % 500 E 1, PONGO t_1000=0.
t=zeros(M,1); % INIZIALIZZAZIONE "t".
for n=M:-1:2
    j=n-1;
    t(j)=(1-t(n))/n;
end

%--------------------------------------------
% ESEGUE UNA MATRICE 2 x 1 DI GRAFICI.
% IL PRIMO GRAFICO LO METTE IN (1,1).
%--------------------------------------------
subplot(2,1,1); 

%------------------------------------------
% PRIMO GRAFICO: PLOT SEMI-LOGARITMICO. 
%------------------------------------------
hold on;
semilogy(1:length(s),abs(s),'k-'); 
semilogy(1:length(s),abs(t(1:length(s))),'m-');
hold off;


%--------------------------------------------
% ERRORI ASSOLUTI.
% ANALISI DI t(1) POSTO t(m)=0; m=2,4,...,20
%--------------------------------------------
t_1_exact=exp(-1);
for m=1:10
    
    M(m)=2*m;
    tt=zeros(M(m),1); % INIZIALIZZAZIONE "t".
    
    for n=M(m):-1:2  % SI OSSERVI CHE IL CICLO VA DA "M" A 2.
        tt(n-1)=(1-tt(n))/n;
    end
    
    val_t_1(m)=abs(t_1_exact-tt(1));
    fprintf('\n \t [M]: %2.0f [VAL.]: %10.15f',M(m),val_t_1(m));
    
end

%--------------------------------------
% SECONDO GRAFICO: ERRORI ASSOLUTI.
%--------------------------------------
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). Sono stati calcolati i valori assoluti perch\mye la successione in avanti assume valori negativi, non permettendo cos\myi l'uso di {\tt{semilogy}}}
 \end{figure}

\noindent
Alcune osservazioni sul codice
{\smsp}
\begin{itemize}
\item Nel primo grafico sono stati calcolati i valori assoluti di $s_n$ perch\mye la successione in avanti assume valori negativi, non permettendo cos\myi l'uso di {\tt{semilogy}} nel plot successivo;
\smsp
\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}

\subsection{Analisi dei risultati ottenuti}
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$.

 \begin{figure}\label{fig1s1}
 \centering
   {\includegraphics[scale=0.4,clip]{grafico_indietro_1_20.eps}}
   {\includegraphics[scale=0.4,clip]{grafico_indietro_1_1000.eps}}
 \caption{Grafico che illustra la successione in avanti: i primi 20 e i primi 1000 termini. Si ricordi che $\exp{(-1)} \approx 0.3679$.}
 \end{figure}

\noindent
Il grafico {\ref{fig3}} mostra per $m=2,\ldots,20$ il comportamento del metodo all'indietro nell'appossimare il valore iniziale
$I_1=\exp{(-1)}$. Evidentemente, al crescere di $m$, l'approssimazione diventa sempre pi\myu accurata. In particolare l'errore assoluto compiuto \mye stampato
nel workspace 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 l'errore assoluto compiuto dalla successione
all'indietro nell'approssimare $\exp{(-1)}$ (partendo da $t_2=0,\ldots,t_20=0$).}
 \end{figure}




\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.0131607401295249$$
%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.0131607401295249
%$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 sul condizionamento}
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)|.$$  Di conseguenza \mye importante valutare la quantit\mya
$$ {\cal{K}}(f,x,x_c)=\frac{ |f(x)-f(x_c)|/|f(x)}{|x-x_c|/|x|}.$$
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},$$ ricaviamo facilmente che
$$ {\cal{K}}(f,x,x_c) \approx {\cal{K}}(f,x) := \frac{|  x \cdot f^{\prime}(x)|}{|f(x)|} .$$


Date le funzioni
$$f_1(x)=1-\sqrt{1-x^2}$$
$$f_2(x)=1-x$$
si calcoli analiticamente il {\sl{condizionamento}}
$${\cal{K}}(f_1,x)=\frac{|  x \cdot f_1^{\prime}(x)|}{|f_1(x)|} $$
$${\cal{K}}(f_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}}(f_1,\cdot)$ e ${\cal{K}}(f_2,\cdot)$ 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{Esercizi proposti}
{\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\mya,
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\myo dire che anche
{\sl{numericamente}} $\tan(\arctan(x))=x$?
{\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$.

% ESERCIZIO 8.
\item {\bf{Esercizio, facoltativo}}.
Per il calcolo della seguente successione di integrali definiti
\[ I_{n} = \int_0^1 \frac{x^n}{x+5} \, dx, \]
si consideri la seguente formula ricorsiva per il calcolo di $I_n$:
\begin{eqnarray*}
s_0  &=& \ln (1.2)  \\
s_{n}  &=&  \frac{1}{n} - 5 s_{n-1} \qquad n > 0 
\end{eqnarray*}
Considerando che la successione $I_n$ \`e descrescente e verifica $0 < I_n <
1/6, \ \forall n$,
si scriva un m-file che implementa l'algoritmo precedente e calcola $s_n$
per un $n$ fissato.
Si calcoli ad esempio $I_{40}$. La ricorrenza \`e stabile?
Per calcolare $I_{40}$ utilizzare un algoritmo stabile. \\

\noindent
{Suggerimento}.Per l'algoritmo stabile si parta ad esempio da
un valore approssimato di $I_{65}$.
{\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}
\end{enumerate}


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



\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{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{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{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{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{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}




