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





\begin{document}

\title{\bf Equazioni nonlineari \thanks{Ultima revisione: \today}}

\author{{A. Sommariva
{\thanks{Dipartimento di Matematica Pura ed Applicata, Universit\`a degli Studi di  Padova, via Trieste 63, 35121 Padova, Italia ({\tt alvise@math.unipd.it})}}} {\mbox{ e}}  {M.Venturin
{\thanks{Dipartimento di Informatica, Universit\`a degli Studi di Verona, Strada Le Grazie 15, 37134 Verona, Italia ({\tt manolo.venturin@gmail.com})}} }
}

\maketitle
\maketitle
\section{Soluzione numerica di equazioni nonlineari}

Data una funzione continua $y=f(x)$, si desidera calcolare $x^{\ast}$ tale
che $f(x^{\ast})=0$. Questo problema \`e assai diffuso nel calcolo
numerico e richiede in generale l'utilizzo di un {\sl{metodo iterativo}}
per approssimare tali {\sl{soluzioni}} $x^{\ast}$. In altri termini, partendo da un valore
iniziale $x_{0}$ si
genera una sequenza di valori $x_{1},x_{2},x_{3},\ldots$ che si desidera
convergano, magari velocemente, ad una opportuna soluzione $x^{\ast}$.

Un esempio significativo \`e rappresentato dal caso in
cui $f$ sia un polinomio $p_{n}$ di grado $n$, ed \`e ben noto da un
teorema
di Galois che per $n \geq 5$, non esistono formule risolutive per il
calcolo degli zeri dell'equazione $p_{n}(x)=0$ che richiedano un
numero finito di operazioni. Ci\myo non vieta di approssimare tali radici
con un metodo numerico
compiendo un errore assoluto e/o relativo inferiore a un valore
prestabilito detto {\sl{tolleranza}}.

Osserviamo che quelle polinomiali non sono gli unici esempi di equazioni
nonlineari. Si consideri ad esempio il problema
di calcolare gli zeri di $f(x)=\sin x -x$, cio\mye quei valori $ x^{\ast}
$ per cui $\sin x^{\ast}-x^{\ast}=0$. Visto che
$$\sin x^{\ast}-x^{\ast}=0 \, \Leftrightarrow \, \sin x^{\ast}= x^{\ast} $$
e
$$|\sin x| \leq 1$$
risulta evidente che se esiste un tale zero, necessariamente
$$x^{\ast} \in [-1,1].$$

Eseguiamo il codice Matlab/Octave
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> x=-10:0.01:10;
>> y=sin(x)-x;
>> plot(x,y);

\end{Verbatim}
Dal plot appare chiaro che se esiste uno zero in $[-10,10]$ allora sta
nell'intervallo $[-1,1]$. Scriviamo quindi
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> x=-1:0.01:1;
>> y=sin(x)-x;
>> plot(x,y);

\end{Verbatim}
L'unico zero sembra essere $x^{\ast}=0$. In effetti la funzione $f$ \mye
continua e decrescente essendo $f^{\prime}(x)=\cos x -1 \leq 0$ (poich\mye
$|\cos x| \leq 1$) ed assume valori positivi e negativi. Quindi ha un
unico zero $x^{\ast}=0$ (si osservi che essendo $\sin 0=0$ abbiamo $\sin
0-0=0$).

In un metodo iterativo, ci sono due aspetti specifici di
cui \`e opportuno tenere conto
\smsp
\begin{enumerate}
\item Garanzia della convergenza alla soluzione: se $\{x_{n}\}$ \`e la
soluzione generata dal metodo e $x^{\ast}$ \`e uno zero per $f$, si
cercano delle condizioni per cui $x_{n} \rightarrow x^{\ast}$.\smsp
\item Se $x_{n} \rightarrow x^{\ast}$ si cerca la velocit\mya con cui
ci\myo accade. In tal senso, \`e importante calcolare il cosidetto
{\sl{ordine di convergenza}}. Sia $\{x_{k}\}$ una successione convergente
ad $x^{\ast}$ e sia
$e_{k}=x_k-x^{\ast}$ l'errore al passo $k$. Se esiste un numero
$p>0$ e una costante $C \neq 0$ tale che
$$\lim_{n \rightarrow \infty}{ \frac {|e_{k+1}|} {|e_{k}|^{p}} }=C$$
allora $p$ \`e chiamato {\sl{ordine di convergenza}} della successione e
$C$ \`e la {\sl{costante asintotica di errore}}.
\smsp
Vediamo con un semplice esempio di confrontare in Matlab il caso in cui sia per $k=0,1,2,\ldots$
$${ \frac {|e_{k+1}|} {|e_{k}|^{p}} }={ \frac {1} {10} }, \,\, {|e_{0}|=1 }$$
\noindent
cio\mye
$${ {|e_{k+1}|}  }={ \frac {1} {10} } \, {|e_{k}|^{p}}, \,\, {|e_{0}|=1 }.$$
\smsp
Evidentemente
$$\lim_{n \rightarrow \infty}{ \frac {|e_{k+1}|} {|e_{k}|^{p}} }={ \frac {1} {10} },$$
\noindent e il metodo ha ordine di convergenza $p$.
\smsp
Se abbiamo $p=1$ la successione prodotta sar\mya:
$$e_0=1, \, e_1=1/10, \, e_2=1/100, \, e_3=1/1000, \, \ldots $$
mentre per $p=2$ abbiamo facilmente
$$e_0=1, \, e_1=1/10, \, e_2=1/1000, \, e_3=1/10^7, \ldots,$$
per cui si vede che a parit\mya di $C$, maggiore \mye $p$ allora minore \mye l'errore $e_k$ compiuto a parit\mya di $k$.
\smsp
Per rendercene maggiormente conto digitiamo in Matlab:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

C=1/10;
e1(1)=1;
e2(1)=1;

M=6;
for index=1:M
  e1(index+1)=C*e1(index);
  e2(index+1)=C*(e2(index))^2;
end

L=length(e1);
semilogy(0:L-1,e1,'rd-'); hold on; semilogy(0:L-1,e2,'ko-');

\end{Verbatim}
\smsp
\noindent
ottenendo il grafico in figura.

 \begin{figure}[!h]\label{tangenti}
 \centering
   {\includegraphics[scale=0.5,clip]{velocita.eps}}
 \caption{Grafico che illustra l'errore di un metodo con convergenza $p=1$ (in rosso a rombi) e $p=2$ (in nero a cerchietti), per $C=1/10$ ed $e_0=1$.}
 \end{figure}

Dal punto di vista pratico per avere un'idea dell'ordine di convergenza si
effettua un plot semi-logaritmico dell'errore. Se esso appare (circa) come
una retta decrescente vuol dire che la (possibile) convergenza ha ordine
$1$ (e viene detta lineare). Altrimenti se \mye al di sotto di qualsiasi
retta (con coefficiente angolare negativo) si dice {\sl{superlineare}}, se
al di sopra {\sl{sublineare}}.
\end{enumerate}

\section{Metodo di bisezione}
Sia $f:[a,b] \rightarrow {\mathbb R}$ una funzione continua e
supponiamo $f(a)\cdot f(b)<0$. Di conseguenza per il {\sl{teorema
degli zeri di una funzione continua}} l'intervallo $(a,b)$ contiene
almeno uno zero di $f$. Supponiamo per semplicit\mya che sia uno solo,
come nel caso in cui $f$ sia crescente o decrescente.

Il metodo di bisezione pu\myo essere definito in forma algoritmica (cf.
{\cite{COM}}, p.
408):
\smsp
\begin{enumerate}
\item si genera una successione di intervalli ($a_{k}$,$b_{k}$) con
$$f(a_{k})\cdot f(b_{k})<0, $$
$$[a_{k},b_{k}] \subset [a_{k-1},b_{k-1}],$$
$$|b_{k}-a_{k}|={\frac{1}{2}} |b_{k-1}-a_{k-1}|.$$
\item Date due tolleranze $\epsilon_{1}>0$, $\epsilon_{2}>0$ si arresta
l'algoritmo o quando
$$|b_{k}-a_{k}| \leq \epsilon_{1}$$
o quando
\begin{equation}\label{residuocriterio}
|f( (a_{k}+b_{k})/2 )| \leq \epsilon_{2}
\end{equation}

o quando $k > {\mbox{nmax}}$, ove ${\mbox{nmax}}$ \`e un numero massimo di iterazioni. Successivamente descriveremo ed implementeremo una variante di
(\ref{residuocriterio}) all'interno del nostro criterio di arresto.
\end{enumerate}

Fissata una tolleranza $\epsilon$, per
avere un'errore assoluto sulla soluzione inferiore ad $\epsilon$
necessitano al pi\'u
$$n \geq {\frac {\log(b-a)-\log \epsilon}{\log 2}}$$
iterazioni del metodo.

\subsection{Implementazione Matlab del metodo di bisezione}

Forniamo ora un'implementazione Matlab del metodo di bisezione per
calcolare uno zero $x^{\ast}$ di $f(x)=0$ con $f:[a,b] \rightarrow
\R$ funzione continua.
Tale procedura che chiameremo {\tt{bisezione}} \mye definita come segue
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function [c,k,semiampiezza,wres]=bisezione(a,b,tolintv,tolres,maxit,f)

\end{Verbatim}
Le variabili di output sono
\smsp
\begin{enumerate}
\item {\tt{c}} \mye un vettore relativo alla sequenza di approssimazioni della soluzione
dell'equazione nonlineare $f(x)=0$ fornita dal metodo di bisezione;
\smsp
\item {\tt{k}} l'iterazione in cui termina il metodo di bisezione;
\smsp
\item {\tt{semiampiezza}} \mye un vettore contenente la semi-ampiezza di
ogni
intervallo
$[a_j,b_j]$ analizzato dal metodo di bisezione;
\smsp
\item {\tt{wres}} \mye un vettore la cui componente $j$-sima consiste del
residuo pesato
$$ | f(c_j) \, \cdot \,  ( b_j-a_j)/(
f(b_j)-f(a_j) ) |:=|{{f(c_j)} \cdot  {        w_j              }}|   $$
con $c_j$ punto medio di $[a_{j-1},b_{j-1}]$ (il punto
$c_{j}$
coincide con $a_{j}$ oppure $b_{j}$) e
$$w_j:=\left( {\frac{ f(b_j)-f(a_j)  }{( b_j-a_j) }} \right)^{-1}.$$

Se la funzione \mye ripida, allora $|w_j|$ \mye piccolo e quindi visto che il residuo
$$|f(x_{j})|$$
non \mye un buon test (il residuo pu\myo essere grande pure se siamo vicini alla soluzione), lo pesa tenendo conto della pendenza;
se invece la funzione \mye molto piatta nuovamente il residuo non \mye un buon test (il residuo pu\myo essere piccolo pure se siamo lontani dalla soluzione) nuovamente lo pesa tenendo conto che $|w_j|$ \mye grande;


\end{enumerate}
\smsp
mentre quelle di input
\smsp
\begin{enumerate}
\item {\tt{a}} estremo minore dell'intervallo contenente lo zero di $f$;
\smsp
\item {\tt{b}} estremo maggiore dell'intervallo contenente lo zero di
$f$;
\smsp
\item {\tt{tolintv}} massima semi-ampiezza ammissibile dell'ultimo
intervallo
analizzato dal metodo di bisezione;
\smsp
\item {\tt{tolres}} tolleranza del residuo pesato;
\smsp
\item {\tt{maxit}} numero massimo di iterazioni del metodo di bisezione;
\smsp
\item {\tt{f}} funzione di cui si vuole calcolare l'unico zero contenuto
in $[a,b]$; la funzione $f$ deve essere
registrata in un m-file come ad esempio
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function y=f(x)
y=sin(x)-x;

\end{Verbatim}
\end{enumerate}

\noindent
Un codice che implementa il metodo di bisezione \mye il seguente:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function [c,k,semilunghezza,residuopesato]=bisezione(a,b,tolintv,...
    tolres,maxit,f)

%--------------------------------------------------------------------------
% NOTA 1.
%
% SE L'UTENTE VUOLE PASSARE LA FUNZIONE g SALVATA IN g.m, SCRIVERA'
%    [...]=bisezione(...,maxit,'g').
%
%--------------------------------------------------------------------------
% NOTA 2.
%
% SI PUO' DEFINIRE LA FUNZIONE 'g' ANCHE USANDO "inline".
%
% ESEMPIO:
%
% >> g=inline('sin(x)-x');
% >> [c,k,semilunghezza,residuopesato]=bisezione(-0.45,0.4,10^(-6),...
%    10^(-6),40,'g');
%	 [IT]:  1 [c]: -0.02500 [a]:-0.025[AMP]: 2.13e-01 [WRES]:1.54e-03
%    ...
%	 [IT]: 18 [c]: 0.40000 [a]:0.400[AMP]: 1.62e-06 [WRES]:1.85e-01
% >>
%--------------------------------------------------------------------------


% PONI "a < b" ANCHE IN CASO DI ERRORE DELL'UTENTE.
if b < a
    s=b; b=a; a=s;
end

fa=feval(f,a);
fb=feval(f,b);

% SE UNO TRA a E b E' LA SOLUZIONE, ALLORA ESCI DOPO AVER
% ASSEGNATO L'OUTPUT.

if fa == 0
    c=a; k=0; semilunghezza=(b-a)/2; residuopesato=0;
    return;
end

if fb == 0
    c=b; k=0; semilunghezza=(b-a)/2; residuopesato=0;
    return;
end

% SE IL NUMERO DI ITERAZIONI NON E' ECCESSIVO ALLORA PROCEDI.

for index=1:maxit

    c(index)=(a+b)/2;
    fc=feval(f,c(index));

    % CALCOLA IL NUOVO INTERVALLO [a,b].
    if sign(fc) == sign(fa)
        a=c(index); fa=fc;  %SUB INTV. DX.
    else
        b=c(index); fb=fc;  %SUB INTV. SX.
    end

    % CALCOLA SEMILUNGHEZZA E RESIDUO PESATO.
    semilunghezza(index)=(b-a)/2;
    den=(fb-fa);
    if den == 0
       den=eps; % SE IL DENOMINATORE "den" E' NULLO, PONGO "den=eps".
    end
    w=(b-a)/den;
    residuopesato(index)=abs(fc*w);

    % SE QUALCHE TEST DI ARRESTO E' VERIFICATO, ESCI DALLA FUNZIONE.
    if (residuopesato(index) < tolres) | ...
             (semilunghezza(index) < tolintv) | (fc == 0)
        k=index;
        fprintf('\n');
        return;
    end

   fprintf('\n \t [IT]:%3.0f [c]: %5.5f [a]:%3.3f',index,c(index),a)
   fprintf('[AMP]: %2.2e [WRES]:%2.2e',semilunghezza(index),...
        residuopesato(index));

end

k=maxit;

fprintf('\n')

\end{Verbatim}
Descriviamo passo passo il programma (a meno di dettagli che lasciamo al
lettore).
\smsp
\begin{enumerate}
\item Supposto $ I_0=[a_0,b_0] \equiv [a,b] $ con
$f(a)\,\cdot\,f(b)
< 0$,
calcoliamo $ c_1 $ punto medio di $I_0$. Se $f(c_1)=0$ si esce,
determinando i parametri richiesti dal metodo.
Altrimenti si determina l'intervallo $ I_1=[a_1,b_1] $ con $
a_1 = c_1 $ e $ b_1 = b_0 $ oppure $ b_1 =
c_1 $ e $ a_1 = a_o $ in modo che sia $f(a_1)\, \cdot \,
f(b_1) < 0$.
\smsp
\item Registrata la semi-ampiezza di $I_1$ nella prima componente del
vettore {\tt{semilunghezza}} e il residuo pesato
$$\rho_{1}:= |f(c_1) \cdot \,  ( b_1 - a_1 )/(
f(b_1)-f(a_1) )| $$ nella prima componente del vettore
{\tt{residuopesato}}, si testa se i criteri di arresto sono verificati,
cio\mye se la semi-ampiezza dell'intervallo \mye inferiore a
{\tt{tolintv}} o il residuo pesato \mye minore di {\tt{tolres}} o si
\mye raggiunto un numero massimo di iterazioni {\tt{maxit}}. Se uno solo
di questi test viene soddisfatto si esce dalla routine altrimenti si
continua.
\smsp
\item Supposto $ I_k=[a_{k},b_{k}]$ con
$f(a_k)\,\cdot\,f(b_k) <
0$,
calcoliamo $ c_{k+1} $ punto medio di $I_k$. Se $f(c_{k+1})=0$ si
esce, determinando i parametri richiesti dal metodo.
Altrimenti si considera l'intervallo $ I_{k+1} =[a_{k+1},b_{k+1}]
$
con $ a_{k+1} = c_{k+1} $ e $ b_{k+1} = b_{k} $ oppure $
b_{k+1} = c_{k+1} $ e $ a_{k+1} = a_{k} $ in modo che sia
$f(a_{k+1})\, \cdot \, f(b_{k+1}) < 0$.
\smsp
\item Registrata la semi-ampiezza di $I_{k+1}$ nella $k+1$-sima componente
del vettore {\tt{semiamp}} e il residuo pesato
$$ \rho_{k+1}:=   |f(c_{k+1}) \cdot \,  ( b_{k+1} -a_{k+1} )/(
f(b_{k+1})-f(a_{k+1}) )| $$ nella $k+1$-sima componente del
vettore
{\tt{wres}}, si testa se i criteri di arresto sono verificati,
cio\mye se la semi-ampiezza dell'intervallo $I_{k+1}$ \mye inferiore a
{\tt{tolintv}} o il residuo pesato \mye minore di {\tt{tolresiduo}} o si
\mye raggiunto un numero massimo di iterazioni {\tt{maxit}}. Se uno solo
di questi test viene soddisfatto si esce dalla funzione altrimenti si
continua.
\end{enumerate}
\smsp
Notiamo alcune cose:
\smsp
\begin{enumerate}
\item Per far eseguire ai programmi esattamente $50$ iterazioni, si pu\myo
impostare il numero massimo di iterazioni a $50$ e le tolleranze sulle
semi-ampiezze e i residui uguali a $0$ (o perfino a $-1$).
\smsp
\item si supponga sia {\tt{f}} una funzione
registrata in un m-file come ad esempio
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function y=f(x)
y=sin(x)-x;

\end{Verbatim}
\noindent
e la si voglia valutare in un punto $x$. La chiamata
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

y=feval('f',x)

\end{Verbatim}
\noindent
valuta la funzione $f$ nel punto $x$. Se $f$ \mye funzione vettoriale, allora $x$ pu\myo essere un vettore. Vediamo un esempio, supposto che in $f$ sia registrata la funzione sopra menzionata.
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> y=feval('f',x)

y =

   -3.1416
>> y=feval(f,x)
??? Input argument 'x' is undefined.

Error in ==> C:\MATLAB6p1\work\f.m
On line 2  ==> y=sin(x)-x;
>>

\end{Verbatim}
\noindent
Notiamo che la mancanza di apici in $f$ genera un errore di sintassi. Nella chiamata della procedura di bisezione, la variabile $f$ viene passata  come stringa '$f$' e quindi all'interno della procedura di bisezione, $f$ \mye per l'appunto una stringa, da cui si pu\myo scrivere
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

fa=feval(f,a);
fb=feval(f,b);

\end{Verbatim}
\noindent
invece di
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

fa=feval('f',a);
fb=feval('f',b);

\end{Verbatim}
\smsp
\item A volte pu\myo essere comodo definire una funzione $g$ senza scriverla su un file {\tt{g.m}}, bens\myi utilizzando il comando
{\tt{inline}}. Un esempio \mye il seguente
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> g=inline('sin(x)-x');
>> [c,k,semilunghezza,residuopesato]=bisezione(-0.45,0.4,...
        10^(-6),10^(-6),40,'g');

   [IT]:  1 [c]: -0.02500 [a]:-0.025[AMP]: 2.13e-01 [WRES]:1.54e-03
   [IT]:  2 [c]: 0.18750 [a]:0.188[AMP]: 1.06e-01 [WRES]:5.52e-02
          ... ... ...
   [IT]: 17 [c]: 0.39999 [a]:0.400[AMP]: 3.24e-06 [WRES]:1.85e-01
   [IT]: 18 [c]: 0.40000 [a]:0.400[AMP]: 1.62e-06 [WRES]:1.85e-01
>>

\end{Verbatim}


\smsp
\item Per disegnare il grafico degli errori relativi a {\tt{rad2}}
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function [fx]=rad2(x)
fx=x^2-2;

\end{Verbatim}
\noindent
si esegua
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

[x,k,semiamp,wres]=bisezione(1,2,0,0,50,'rad2');
x_fin=x(length(x));
semilogy(abs(x-x_fin)/x_fin,'r-');

\end{Verbatim}
Tale codice produce il grafico degli errori relativi al metodo. Osserviamo che
{\tt{semilogy}} chiede nell'ordine
il vettore di ascisse e il vettore di ordinate. Qualora ci sia un solo
vettore lo prende quale vettore delle ordinate, con ascisse il
vettore degli indici da $1$ alla lunghezza del vettore delle ordinate.
\end{enumerate}
\smsp
\noindent
{\bf{Da provare}} Salvare su un file {\bf{f.m}} la function
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function fx=f(x)
fx=sin(x)-x;

\end{Verbatim}
\noindent
e quindi digitare dal workspace di Matlab/Octave
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

[x,k,semiamp,wres]=bisezione(-3,2,0,0,50,'f');

\end{Verbatim}
Il risultato come previsto \mye $x=0$. Pi\myu precisamente
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

 [IT]: 1 [c]: -0.50000 [z1]:-0.500 [AMP]:1.25e+000 [WRES]:4.63e-002
 [IT]: 2 [c]:  0.75000 [z1]:-0.500 [AMP]:6.25e-001 [WRES]:9.61e-001
 [IT]: 3 [c]:  0.12500 [z1]:-0.500 [AMP]:3.13e-001 [WRES]:9.73e-003
 [IT]: 4 [c]: -0.18750 [z1]:-0.188 [AMP]:1.56e-001 [WRES]:2.41e-001
 [IT]: 5 [c]: -0.03125 [z1]:-0.031 [AMP]:7.81e-002 [WRES]:2.41e-003
 [IT]: 6 [c]:  0.04688 [z1]:-0.031 [AMP]:3.91e-002 [WRES]:6.03e-002
 [IT]: 7 [c]:  0.00781 [z1]:-0.031 [AMP]:1.95e-002 [WRES]:6.01e-004
 [IT]: 8 [c]: -0.01172 [z1]:-0.012 [AMP]:9.77e-003 [WRES]:1.51e-002
 [IT]: 9 [c]: -0.00195 [z1]:-0.002 [AMP]:4.88e-003 [WRES]:1.50e-004
 [IT]:10 [c]:  0.00293 [z1]:-0.002 [AMP]:2.44e-003 [WRES]:3.77e-003
 ......
 [IT]:40 [c]: -0.00000 [z1]:-0.000 [AMP]:2.27e-012 [WRES]:4.55e-012
 [IT]:41 [c]: -0.00000 [z1]:-0.000 [AMP]:1.14e-012 [WRES]:2.27e-012
 [IT]:42 [c]: -0.00000 [z1]:-0.000 [AMP]:5.68e-013 [WRES]:0.00e+000
 [IT]:43 [c]: -0.00000 [z1]:-0.000 [AMP]:2.84e-013 [WRES]:5.68e-013
 [IT]:44 [c]: -0.00000 [z1]:-0.000 [AMP]:1.42e-013 [WRES]:0.00e+000
 [IT]:45 [c]: -0.00000 [z1]:-0.000 [AMP]:7.11e-014 [WRES]:0.00e+000
 [IT]:46 [c]: -0.00000 [z1]:-0.000 [AMP]:3.55e-014 [WRES]:0.00e+000
 [IT]:47 [c]: -0.00000 [z1]:-0.000 [AMP]:1.78e-014 [WRES]:3.55e-014
 [IT]:48 [c]: -0.00000 [z1]:-0.000 [AMP]:8.88e-015 [WRES]:0.00e+000
 [IT]:49 [c]: -0.00000 [z1]:-0.000 [AMP]:4.44e-015 [WRES]:8.88e-015
 [IT]:50 [c]: -0.00000 [z1]:-0.000 [AMP]:2.22e-015 [WRES]:4.44e-015

\end{Verbatim}
\noindent
Notiamo che a volte sul workspace compare
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

error: `c' undefined near line 3 column 8
error: evaluating argument list element number 1
error: evaluating binary operator `-' near line 3, column 10
error: evaluating assignment expression near line 3, column 3
error: ...

\end{Verbatim}
\noindent
In molti casi \mye sufficiente rilanciare il comando
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

[c,k,semiamp,wres]=bisezione(-3,2,0,0,50,'f')

\end{Verbatim}
\noindent
per avere i risultati desiderati.

\section{Esercizio per casa} Calcolare con i metodi di bisezione e Newton un'approssimazione dello zero $\alpha$ di
\begin{eqnarray}
&& 5x = {\exp{(-x)}}  \\
&& x^4-2 = 0 \\
&& x \log{(x)} -1 =0\\
&& (x-1) \log{(x)} =0 \\
&& m  = x - E {\sin{(x)}}, \,\, m=0.8, \, E=2
\end{eqnarray}
Verificare che per le approssimazioni $x^{\ast} \approx\alpha$ fornite sia effettivamente $f({x^{\ast}}) \approx 0$. Posto come $x^{\ast}$ il risultato finale, si disegni il grafico degli errori $|x^{(k)}-x^{\ast}|$ in scala semilogaritmica. Da questa dedurre, se possibile, l'ordine di convergenza del metodo e giustificare teoricamente la risposta.
\smsp
Alcuni suggerimenti:
\smsp
\begin{enumerate}
\item si controlli sempre che la directory corrente contenga i files Matlab/Octave da usare; ad esempio se voglio usare {\tt{bisezione}} e una funzione {\tt{f.m}} devo controllare che siano presenti nella directory corrente come nell'esempio
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> ls
Dg1.m   dg.m   g.m   rad.m
bisezione.m   drad2.m   g1.m   rad2.m
bisezione.m~   f.m   newton.m
>>

\end{Verbatim}
A volte appaiono files con la {\sl{tilda}}, che sono alcuni backups automatici di Matlab.
\smsp
\item
alcuni utenti di Octave possono incontrare problemi con l'editor {\tt{Kate}} che aggiunge caratteri indesiderati al testo; quale soluzione si utilizzi invece {\tt{pico}}di {\tt{Kate}};
\smsp
\item si ricordi che la chiamata di bisezione \mye del tipo
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

[c,k,semiamp,wres]=bisezione(-3,2,0,0,50,'f')

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

[c,k,semiamp,wres]=bisezione(-3,2,0,0,50,f)

\end{Verbatim}
(si scriva cio\mye {\tt{'f'}} e non {\tt{f}} per passare la funzione da studiare).
\smsp
\item pu\myo essere utile descrivere la funzione {\tt{f}} come inline; ad esempio
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> f=inline('sin(x)-x');
>> [c,k,semiamp,wres]=bisezione(-0,35,0.45,0.001,0.001,50,f)

\end{Verbatim}
\smsp
\item si osservi che alcune funzioni non sono scritte come $f(x)=0$ che \mye la formulazione richiesta; a tale scopo portare tutti i termini a primo membro e determinare la $f$.
\smsp
\item per determinare gli estremi di bisezione $a$, $b$ a volte pu\myo essere utile eseguire il grafico della funzione.
\end{enumerate}

\smsp

\begin{remark}
Il metodo di bisezione non ha, secondo le definizioni date, ordine di convergenza 1, come molti studenti ritengono. Vediamo il seguente test, relativamente all'equazione
$$\sin{(x)}-x=0.$$
Apriamo il file {\tt{test\_bisezione}}
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]


function ratio=test_bisezione(esempio)

% BISEZIONE E ORDINE DI CONVERGENZA.

switch esempio
case 1
    f=inline('sin(x)-x'); % HA ZERO IN "sol=0".
    sol=0;
    a=-3; b=2;
case 2
    f=inline('x.^2-2'); % HA ZERO IN "sqrt(2)".
    sol=sqrt(2);
    a=1; b=2;
end

tolintv=10^(-6);
tolres=10^(-6);
maxit=20;
[c,k,semilunghezza,residuopesato]=bisezione(a,b,tolintv,...
    tolres,maxit,f);

en=abs(sol-c);
e_n=en(1:end-1); % e_n=abs(sol-c_n)
e_nplus1=en(2:end);   % e_nplus1=(sol-c_n+1)

% RIDUZIONE n-SIMA ITERAZIONE.
ratio=abs(e_nplus1./e_n);

\end{Verbatim}
Studiamo il primo caso, cio\mye quello relativo al problema $\sin{(x)}-x=0$.
Per prima cosa osserviamo che $0$ \mye la soluzione dell'equazione $\sin{(x)}-x=0$ \mye $x^{\ast}=0$. Prendiamo quale intervallo iniziale $[a,b]=[-3,2]$.
Lanciamo il programma appena descritto e otteniamo
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> ratio=test_bisezione(1)

   [IT]:  1 [c]: -0.50000 [a]:-0.500[AMP]: 1.25e+00 [WRES]:4.63e-02
   [IT]:  2 [c]: 0.75000 [a]:-0.500[AMP]: 6.25e-01 [WRES]:9.61e-01
   [IT]:  3 [c]: 0.12500 [a]:-0.500[AMP]: 3.12e-01 [WRES]:9.73e-03
   [IT]:  4 [c]: -0.18750 [a]:-0.188[AMP]: 1.56e-01 [WRES]:2.41e-01
   [IT]:  5 [c]: -0.03125 [a]:-0.031[AMP]: 7.81e-02 [WRES]:2.41e-03
   [IT]:  6 [c]: 0.04688 [a]:-0.031[AMP]: 3.91e-02 [WRES]:6.03e-02
   [IT]:  7 [c]: 0.00781 [a]:-0.031[AMP]: 1.95e-02 [WRES]:6.01e-04
   [IT]:  8 [c]: -0.01172 [a]:-0.012[AMP]: 9.77e-03 [WRES]:1.51e-02
   [IT]:  9 [c]: -0.00195 [a]:-0.002[AMP]: 4.88e-03 [WRES]:1.50e-04
   [IT]: 10 [c]: 0.00293 [a]:-0.002[AMP]: 2.44e-03 [WRES]:3.77e-03
   [IT]: 11 [c]: 0.00049 [a]:-0.002[AMP]: 1.22e-03 [WRES]:3.76e-05
   [IT]: 12 [c]: -0.00073 [a]:-0.001[AMP]: 6.10e-04 [WRES]:9.42e-04
   [IT]: 13 [c]: -0.00012 [a]:-0.000[AMP]: 3.05e-04 [WRES]:9.39e-06
   [IT]: 14 [c]: 0.00018 [a]:-0.000[AMP]: 1.53e-04 [WRES]:2.35e-04
   [IT]: 15 [c]: 0.00003 [a]:-0.000[AMP]: 7.63e-05 [WRES]:2.35e-06
   [IT]: 16 [c]: -0.00005 [a]:-0.000[AMP]: 3.81e-05 [WRES]:5.89e-05
ratio =

  Columns 1 through 14

    1.5000    0.1667    1.5000    0.1667    ....

  Columns 15 through 16

    1.5000    0.1667

>> 

\end{Verbatim}

\noindent La successione di rapporti $|e_{n+1}/e_n|$ ovviamente non ha limite per $n \rightarrow +\infty$ poich\mye assume alternativamente valori $1.5$ e $1/6 \approx 0.1667$. Quindi non si pu\myo affermare che l'ordine di bisezione sia $1$.

Tenendo in mente quanto appena visto, detto $|e_{n+1}|/|e_n|$ il fattore di riduzione alla $n$-sima iterazione, e indicato come fattore di riduzione medio dopo $n$ iterazioni la quantit\mya
$$k_n(e_0)=( {\frac{|e_{n}|}{|e_{n-1}|}} \ldots  {\frac{|e_{1}|}{|e_0|}}   )^{1/n}=( {\frac{|e_{n}|}{|e_{0}|}})^{1/n}$$
verifichiamo che per $n \rightarrow \infty$ il fattore di riduzione medio tende a $1/2$ (il che ovviamente non vuol dire che la riduzione sia uguale al massimo $1/2$ come erroneamente ritenuto da molti studenti).
Per convincersene analiticamente, osserviamo che nel nostro esempio
\smsp
\begin{itemize}
\item se $n$ \mye pari allora per un certo $m$ si ha $n=2m$ ed \mye $$k_{2m}(e_0)=((3/2)^m (1/6)^m)^{1/(2m)}=((3/12)^m)^{1/(2m)}=1/2$$
\smsp
\item se $n$ \mye dispari allora per un certo $m$ si ha $n=2m+1$ ed \mye dopo facili conti
$$k_{2m+1}= (3/2)^{1/(2m+1)} (1/4^m)^{1/(2m+1)} = (3/2)^{1/(2m+1)} 2^{-(2m)/(2m+1)} \rightarrow 1/2.$$
\end{itemize}

\smsp
E' d'altra parte evidentemente che il fattore di riduzione medio nulla ha a che fare con la classica media dei fattori di riduzione. Calcoliamo le medie al variare di $n$ e vediamo che il limite vale $5/6$ e non $1/2$.
\smsp
\begin{itemize}
\item se $n$ \mye pari allora per un certo $m$ si ha $n=2m$ ed \mye
$${\mbox{media}}_{2m}= \frac{m \cdot(3/2)+m \cdot (1/6)}{2m} = \frac{(3/2)+ (1/6)}{2}= 10/12 \approx 0.8334.$$

\smsp
\item se $n$ \mye dispari allora per un certo $m$ si ha $n=2m+1$ ed \mye dopo facili conti
\begin{eqnarray}
{\mbox{media}}_{2m+1} &=& \frac{(m+1) \cdot(3/2)+m \cdot (1/6)}{2m+1} \nonumber \\
&=& \frac{(m+1) \cdot(3/2)}{2m+1} + \frac{m \cdot (1/6)}{2m+1} \rightarrow (1/2)\cdot(3/2)+(1/2) \cdot (1/6) \nonumber \\
&=& 10/12 \approx 0.8334.
\end{eqnarray}
\end{itemize}

\end{remark}



\section{Metodo di Newton per la risoluzione di equazioni nonlineari}

Supponiamo che $f$ sia derivabile con continuit\mya su un sottinsieme di
${\mathbb{R}}$ e sia $f^{(1)}(x)$ la derivata prima di $f$
valutata nel generico punto $x$.

Osserviamo che  se $x^{\ast}$ \mye uno zero di $f$ in $[a,b]$ e se $f$ \mye sufficientemente regolare in $[a,b]$  allora dalla formula di Taylor centrata in $x_k \in [a,b]$, per un certo $\xi_k$ che sta nel pi\myu piccolo intervallo aperto contentente $x_k$, $x^{\ast}$, abbiamo
$$0=f(x^{\ast})=f(x_k)+f^{\prime}(x_k)(x^{\ast}-x_k)+f^{\prime \prime}(\xi_k)(x^{\ast}-x_k)^2/2$$
Di conseguenza, tralasciando i termini di ordine superiore
$$0 \approx f(x_k)+f^{\prime}(x_k)(x^{\ast}-x_k)$$
\noindent 
e quindi se $f^{(1)}(x_{k}) \neq 0$, dopo facili conti
$$x^{\ast} \approx x_{k}-f(x_{k})/f^{(1)}(x_{k}).$$

\smsp

Il {\sl{metodo di Newton}} genera la successione
\begin{equation}\label{newton}
x_{k+1}=x_{k}+h_{k}, \, \, h_{k}=-f(x_{k})/f^{(1)}(x_{k}), \,\,
k=0,1,\ldots
\end{equation}
supposto che sia $f^{(1)}(x_{k}) \neq 0$ per $k=0,1,\ldots$.

\subsection{Alcuni risultati sulla convergenza locale e globale di Newton}
Per quanto riguarda la velocit\mya di convergenza si dimostra il
seguente teorema di {\sl{convergenza locale}} (cf. {\cite{ATK}}, p. 60)
\smsp

\begin{theorem}\label{conv1}
Sia $x^{\ast}\in (a,b)$ uno zero semplice di $f: [a,b] \rightarrow \R$. Si supponga inoltre $f \in C^{2}([a,b])$. Allora per $x_0 \in [a,b]$ sufficientemente vicino a $x^{\ast}$ le iterazioni del metodo di Newton
$$x_{k+1}=x_k-f(x_k)/f^{(1)}(x_k), \,\,\, k=0,1,2,\ldots$$
sono ben definite e convergono quadraticamente a $x^{\ast}$.
\end{theorem}

\smsp

\begin{proof}
La dimostrazione consiste dapprima nel trovare un intervallo ${\cal{I}}_1$ contenente lo zero in cui la derivata esiste e poi un sottointervallo ${\cal{I}}_2$ per cui se $x_n \in {\cal{I}}_2$ allora l'iterazione $x_{n+1}$ del metodo di Newton fa ancora parte di ${\cal{I}}_2$. Da stime dell'errore in quest'intervallo si deduce la convergenza quadratica. Vediamo i dettagli della dimostrazione.

\smsp

Essendo lo zero semplice, per definizione $f^{(1)}(x^{\ast}) \neq 0$. Di conseguenza esiste un intervallo simmetrico chiuso ${\cal{I}}_1=[x^{\ast}-\delta,x^{\ast}+\delta]$ di $x^{\ast}$ in cui $f^{(1)} \neq 0$. Osserviamo che dalla formula di Taylor centrata in $x_n$ e dalla definizione della successione del metodo di Newton, se $x_n \in {\cal{I}}_1$, allora per un certo $\xi_n \in {\cal{I}}_1$ (che in effetti sta nel pi\myu piccolo intervallo aperto contenente $x^{\ast}$ e $x_n$)
$$0=f(x^{\ast})=f(x_n)+f^{(1)}(x_n) (x^{\ast}-x_n)+f^{(2)}(\xi_n) (x^{\ast}-x_n)^2/2$$
e quindi dividendo ambo i membri per $f^{(1)}(x_n) \neq 0$ (poich\mye $x_n \in {\cal{I}}_1$), ricordando la definizione del metodo di Newton
\begin{eqnarray}\label{loc1}
0 &=& {\frac{f(x_n)}{f^{(1)}(x_n)}} + (x^{\ast}-x_n)+{\frac{f^{(2)}(\xi_n)}{2 \cdot f^{(1)}(x_n)}}     (x^{\ast}-x_n)^2 \nonumber \\
&=& x^{\ast} - x_{n+1} + {\frac{f^{(2)}(\xi_n)}{2 \cdot f^{(1)}(x_n)}}     (x^{\ast}-x_n)^2.
\end{eqnarray}
Di conseguenza, posto $e_n=x_n-x^{\ast}$, da ({\ref{loc1}}) abbiamo
$$|e_{n+1}|=| {\frac{f^{(2)}(\xi_n)}{2 \cdot f^{(1)}(x_n)}} |    |e_n|^2.$$
Per il teorema di Weierstrass esiste finito
$$C={ \max_{x \in {\cal{I}}_1} {|f^{(2)}(x)|} } / { \min_{y \in I_1} {|f^{(1)}(y)|} } .$$
Possiamo quindi affermare che se $x_n \in I_1$ allora
$$|e_{n+1}| \leq C |e_n|^2.$$
Supponiamo ora che per $\epsilon >0$ sia $\gamma=\min{(\delta,1/C)}-\epsilon>0$ e che $x_0 \in {\cal{I}}_2 \equiv [x^{\ast}-\gamma,x^{\ast}+\gamma] \subseteq {\cal{I}}_1$ (si osservi che $\gamma < \delta$). Allora essendo $0 \leq C \gamma = C(\min{(\delta,1/C)}-\epsilon) \leq 1 -C\epsilon < 1$ ricaviamo
$$|e_1| \leq C |e_0|^2 =(C|e_0|) \cdot |e_0| \leq (C \gamma) \gamma < \gamma$$
per cui tutte le iterazioni $x_n$ del metodo di Newton stanno nell'intervallo ${\cal{I}}_2=[x^{\ast}-\gamma,x^{\ast}+\gamma]$ ed essendo ${\cal{I}}_2 \subseteq {\cal{I}}_1$ la successione descritta dal metodo di Newton \mye ben definita.
Anzi, da $x_0 \in {\cal {I}}_2$ implica $|e_0| \leq \gamma$, ed essendo $$0 \leq C|e_{n+1}| \leq (C^2 |e_n|^2)= (C |e_n|)^2 \leq \ldots \leq (C |e_0|)^{2^{n+1}}$$ con $C|e_0| \leq C \gamma <1$,
necessariamente, per il teorema del confronto tra limiti, $|e_{n+1}| \rightarrow 0$ cio\mye $x_{n+1} \rightarrow x^{\ast}$.
Inoltre si ha o che per un certo $n^{\ast}$, se $n \geq n^{\ast}$, \mye $x_n=x^{\ast}$ oppure dalla continuit\mya di $f^{(1)}$ e $f^{(2)}$ nell'intervallo ${\cal{I}}_2=[x^{\ast}-\gamma,x^{\ast}+\gamma]$
$$\lim_n {\frac {|e_{n+1}|}{|e_n|^2}}=\lim_n | {\frac{f^{(2)}(\xi_n)}{2 \cdot f^{(1)}(x_n)}} |={\frac{|f^{(2)}(x^{\ast})|}{2 \cdot |f^{(1)}(x^{\ast})|}} \neq 0   $$
\noindent cio\mye la convergenza \mye quadratica.

\end{proof}

\smsp
Particolarmente interessante il seguente teorema relativo alla {\sl{convergenza}} del metodo di Newton (cf. {\cite[p.62]{ATK}}):
\smsp

\begin{theorem}\label{convatk}
Sia $f \in C^{2}([a,b])$, con $[a,b]$ intervallo chiuso e limitato. Se
\smsp
\begin{enumerate}
\item $f(a) \cdot f(b) < 0$;\smsp
\item $f^{(1)}(x) > 0$, per ogni $x \in [a,b]$;\smsp
\item $f^{(2)}(x) > 0$ \smsp
\end{enumerate}
allora le iterate $x_{k}$ fornite dal metodo di Newton sono strettamente decrescenti e convergono all'unica soluzione $x^{\ast}$ in $[a,b]$, per $x_{0}=b$.
\end{theorem}
\smsp
\begin{remark}
Nelle ipotesi del Teorema \ref{convatk}, dal Teorema \ref{conv1}, si deduce che la convergenza \mye quadratica. Si sottolinea che se queste non sono verificate, non si pu\myo asserire che la convergenza sia quadratica (ad esempio se uno zero $x^{ast}$ \mye multiplo allora la convergenza di Newton \mye solo lineare!).
\end{remark}
\smsp


\begin{remark}
In questa nota consideriamo alcuni casi in cui nonostante il teorema \ref{convatk} non sia direttamente applicabile si possa comunque giungere a provare la convergenza del metodo di Newton.

\smsp
Si supponga sia $f^{(1)}(x) < 0$, $f^{(2)}(x) > 0$ per ogni $x \in [a,b]$ e sia $g(x):=f(-x)$. E' chiaro che se $f(x^{\ast})=0$ si pone i $\tilde{x}=-x^{\ast}$ allora $g({\tilde{x}})=f(-{\tilde{x}})=f(x^{\ast})=0$. Inoltre per $x \in [-b,-a]$ abbiamo $-x \in [a,b]$ e
$$g^{(1)}(x)=-f^{(1)}(-x)>0, \,\, g^{(2)}(x)=f^{(2)}(-x)>0.$$ \noindent Quindi, posto $x^{(0)}=-a$, le ipotesi del teorema valgono per $g(x)$ ed il metodo di Newton produce quindi una successione $\{x_k{}\}_k$ decrescente e convergente all'unica soluzione $\tilde{x}$ di $g(x)=0$ in $[-b,-a]$. Essendo
$$x_{k+1}=x_{k}-g(x_{k})/g^{(1)}(x_{k})=x_{k}-f(-x_{k})/(-f^{(1)}(-x_{k}))=x_{k}+f(-x_{k})/(f^{(1)}(-x_{k}))$$
\noindent
posto $t^{(k)}=-x^{(k)}$ e moltiplicando entrambi i membri della precedente equazione per $-1$, ricaviamo che la successione
$$t_{k+1}=t_{k}-f(t_{k})/(f^{(1)}(t_{k})), \, t_0=a$$
\noindent \mye crescente e convergente all'unica soluzione ${x^{\ast}}$ di $f(x)=0$ in $[a,b]$.

\smsp
Analogamente, se $f^{(1)}(x) < 0$, $f^{(2)}(x) < 0$, al posto di $f(x)=0$ basta studiare l'equazione $-f(x)=0$. Se invece $f^{(1)}(x) > 0$, $f^{(2)}(x) < 0$ per ogni $x \in [a,b]$, allora bisogna studiare invece di $f(x)=0$ l'equazione $h(x):=-f(-x)=0$. Nuovamente se $h({\tilde{x}})=0$ allora $-f(-{\tilde{x}})=0$ da cui $x^{\ast}=-{\tilde{x}}$.
\end{remark}
\smsp


\noindent Dimostriamo una interessante variante del precedente

\smsp

\begin{theorem}
Sia $f: [a,b] \rightarrow \R$ una funzione che gode delle seguente propriet\mya
\begin{enumerate}
\item $f(a)<0$ e $f(b)>0$;
\item $f^{\prime \prime}(x) > 0$ per ogni $x \in [a,b]$.
\end{enumerate}
Allora
\begin{enumerate}
\item esiste uno ed un solo $x^{\ast} \in (a,b)$ tale che $f(x^{\ast})=0$. Tale zero \mye semplice.
\item per $x_0 \in [x^{\ast},b]$ il metodo di Newton
$$x_{n+1}=x_n- {\frac{f(x_n)}{f^{\prime} (x_n) }}$$
\mye ben definito.
\item La successione $\{x_n\}$ \mye decrescente e converge a $x^{\ast}$.
\end{enumerate}
\end{theorem}

\smsp

\begin{proof}

\begin{enumerate}

\item Essendo $f^{\prime \prime}(x) > 0$ per ogni $x \in [a,b]$, allora la funzione $f^{\prime}$ \mye continua e monotona crescente in $[a,b]$. Di conseguenza,
\begin{itemize}
\item se $f^{\prime} > 0$ in $[a,b]$ si ha che pure $f$ \mye monotona crescente§ e quindi in $(a,b)$ c'\mye un unico zero $x^{\ast}$ poich\mye $f(a) <0$ e $f(b)>0$;
\item se $f^{\prime}$ non ha segno costante in $[a,b]$ essendo monotona crescente ha un unico zero. Di conseguenza $f$ ha un unico minimo $x_{\mbox{\tiny{min}}} \in [a,b]$. Quindi dalla monotonia crescente di $f^{\prime}$ in $[x_{\mbox{\tiny{min}}},b]$ e dal fatto che $f(x_{\mbox{\tiny{min}}} ) \leq f(a)$ e $f(x_{\mbox{\tiny{min}}} )  \leq 0 < f(b)$ abbiamo che vi \mye un unico zero in $x^{\ast} \in [a,b]$.
\end{itemize}
Osserviamo che per quanto appena visto, $f^{\prime}(x) > 0$ per ogni $x \in [x^{\ast},b]$ e quindi in particolare che lo zero $x^{\ast}$ \mye semplice.

\item Questo risultato segue dal fatto che la funzione $f^{\prime}$ \mye strettamente positiva in $[x^{\ast},b]$.

\item Sia $n \geq 0$ e $x_n \geq x^{\ast}$. Dalla formula di Taylor, essendo $f^{\prime \prime}(x_n) > 0$ deduciamo che per $\xi_n$ nel pi\myu piccolo intervallo aperto contenente
$x_n$ e $x_{n+1}=x_n- {\frac{f(x_n)}{f^{\prime} (x_n) }}$ si ha
\begin{eqnarray}
f(x_{n+1})&=&f(x_n)+f^{\prime}(x_n) \cdot (x_{n+1}-x_n) +f^{\prime \prime}(\xi_n) \cdot {\frac{(\xi-x_n)^2}{2}} \nonumber \\
&>& f(x_n)+f^{\prime}(x_n) \cdot (x_{n+1}-x_n)=0=f(x^{\ast}).
\end{eqnarray}
Poich\mye per $x < x^{\ast}$ la funzione $f$ \mye negativa, necessariamente ci\myo implica che $x_{n+1} > x^{\ast}$
Inoltre la successione $\{x_n\}$ \mye strettamente decrescente perch\mye da $f(x_n) > 0$ e $f^{\prime} (x_n) >0$ abbiamo $$x_{n+1}=x_n- {\frac{f(x_n)}{f^{\prime} (x_n) }} \leq x_n.$$
Poich\mye la successione $\{x_n\}$ \mye decrescente e limitata inferiormente da $x^{\ast}$, essa ammette limite $\xi^{\ast}$. Ma allora per la continuit\mya di $f$ e $f^{\prime}$
$$\xi^{\ast}=
 \lim_n {x_{n+1}}=  \lim_n {x_n}-\lim_n  {\frac {f(x_n)} {f^{\prime} (x_n)} }= \xi^{\ast} - { \frac {f(\xi^{\ast})}   {f^{\prime} (\xi^{\ast})}      }$$
ed essendo $f^{\prime} (\xi^{\ast}) \neq 0$ poich\mye tale zero \mye semplice abbiamo $f(\xi^{\ast})=0$. Ma $x^{\ast}$ \mye l'unico zero di $f$ e quindi $\xi^{\ast}=x^{\ast}$ cio\mye il metodo di Newton converge decrescentemente a $x^{\ast}$.
\end{enumerate}

\end{proof}

\smsp

\noindent
Una versione modificata di questo teorema che stabilisce un teorema di convergenza globale si trova in {\cite{COM}}, p. 419:
\smsp
\begin{theorem}\label{convglob}
Sia $f \in C^{2}([a,b])$, con $[a,b]$ intervallo chiuso e limitato. Se
\smsp
\begin{enumerate}
\item $f(a) \cdot f(b) < 0$; \smsp
\item $f^{(1)}(x) \neq 0$, per ogni $x \in [a,b]$; \smsp
\item $f^{(2)}(x) \geq 0$ o $f^{(2)}(x) \leq 0$, per ogni $x \in [a,b]$; \smsp
\item $|f(a)/f^{(1)}(a)|<b-a$ e $|f(b)/f^{(1)}(b)|<b-a$, \smsp
\end{enumerate}
allora il metodo di Newton converge all'unica soluzione $x^{\ast}$
in $[a,b]$, per ogni $x_{0} \in [a,b]$.
\end{theorem}
\smsp
Osserviamo che
\smsp
\begin{enumerate}
\item il metodo di Newton non ha sempre convergenza
quadratica, come nel caso del problema $f(x)=0$ in cui $x^{\ast}$ sia uno zero avente molteplicit\mya $p > 1$ cio\mye tale che
$$f(x^{\ast})=f^{(1)}(x^{\ast})=\ldots=f^{(p-1)}(x^{\ast})=0;$$ \noindent quale esempio si ponga $f(x)={\log{(x-1)}}(x-2)$; si verifica facilmente che $f(2)=f^{\prime}(2)=0$ e che la convergenza non \mye quadratica; per convincersene si applichi il metodo di Newton a tale esempio e si esegua il grafico dell'errore compiuto in scala semilogaritmica;
\smsp
\item i due teoremi sopracitati forniscono delle condizioni necessarie per la convergenza; in realt\mya in letteratura ne esistono molte altre (cf. \cite{DSC}, \cite{KEL}, \cite{ZEI}).
\end{enumerate}



\subsection{Implementazione del metodo di Newton}



Dal punto di vista algoritmico quanto visto pu`\myo essere brevemente riassunto dal seguente
pseudo-codice

{\vspace{0.5cm}}

{\textsc{ALGORITMO }} $y=$ {\tt{ newton }}($x$,$\tau$)
\smsp
\begin{enumerate}
\item $y=x$. \smsp
\item Calcola $f^{(1)}(y)$. \smsp
\item Se $f^{(1)}(y)$ non \`e $0$, poni $s =-f(y)/f^{(1)}(y)$. \smsp
\item $y=y+s$.\smsp

\item Do while $\|s\|>\tau$. \smsp
      \begin{enumerate}
       \item Calcola $f^{(1)}(y)$. \smsp
       \item Se $f^{(1)}(y)$ non \`e $0$, poni $s =-f(y)/f^{(1)}(y)$. \smsp
       \item $y=y+s$. \smsp
      \end{enumerate}
\end{enumerate}
dove $x$ \`e il punto iniziale e $\tau$ una tolleranza prefissata. 

Definiamo lo {\sl{scarto}} (assoluto) di un metodo iterativo come
\begin{equation}\label{sc1}
s_{n+1}=x_{n+1}-x_{n}.
\end{equation}
Un test comunemente utilizzato \mye che sia
\begin{equation}\label{sc2}
|s_{n+1}| \leq \tau
\end{equation}
 dove
$\tau$ \mye una tolleranza prefissata dall'utente (ad esempio
$\tau=10^{-6}$).

Implementiamo il metodo di Newton mediante una function, usando come
criteri di arresto lo scarto ({\ref{sc1}})-({\ref{sc2}}) ed il numero massimo di iterazioni. In
seguito testiamo il metodo con diverse combinazioni del punto iniziale
$x_0$, della tolleranza richiesta e del numero massimo di iterazioni e si
commentiamo i risultati. Il codice fornisce, come output, i valori della
successione e dello scarto ad ogni passo nonch\mye il numero di
iterazioni.

\smsp Supponiamo di voler studiare l'equazione $f(x)=0$ con
$$f(x)=\sin{(x)}-x$$
e quindi
$$f^{\prime}(x)=\cos{(x)}-1$$
Scriviamo in files {\tt{f.m}}, {\tt{f1.m}}, {\tt{newton.m}} le
seguenti funzioni:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function y=f(x)
y=sin(x)-x;

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

function y=f1(x)
y=cos(x)-1;

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

function [x,k,scarto,errflag]=newton(x0,tol,kmax,funct,dfunct)

%---------------------------------
% METODO DI NEWTON.
%---------------------------------
% INPUT.
%---------
% x0: punto iniziale.
% tol: tolleranza criterio dello step/residuo pesato.
% kmax: numero massimo di iterazioni.
% funct: funzione di cui cercare lo zero "funct(x)=0".
% funct: derivata di "funct".
%---------------------------------
% OUTPUT.
%---------
% x: vettore iterazioni del metodo di Newton.
% k: ultima iterazione.
% scarto: vettore degli scarti -fx/dfx.
% errflag: 0: uscita corretta, 1: uscita per errore.
%---------------------------------
% ESEMPIO.
%---------
% >>f=inline('x.^2-2');
% >>df=inline('2*x');
% >>[x,k,scarto,errflag]=newton(0.8,10^(-6),40,f,df)
% [ITER.]:   1[VALORE]: 0.80000 [ABS.SCARTO]: 8.50e-01
% [ITER.]:   2[VALORE]: 1.65000 [ABS.SCARTO]: 2.19e-01
% [ITER.]:   3[VALORE]: 1.43106 [ABS.SCARTO]: 1.67e-02
% [ITER.]:   4[VALORE]: 1.41431 [ABS.SCARTO]: 9.92e-05
% [ITER.]:   5[VALORE]: 1.41421 [ABS.SCARTO]: 3.48e-09
% x =
%
%     0.8000
%     1.6500
%     1.4311
%     1.4143
%     1.4142
%
% k =
%
%      5
%
% scarto =
%
%     0.8500
%    -0.2189
%    -0.0167
%    -0.0001
%    -0.0000
%
%  errflag =
%
%      0
%
%---------------------------------
% NOTA.
%---------
% IN ALCUNI MATLAB, SE QUANTO PROVATO DA' ERRORE
% ESEGUIRE
%
% [x,k,scarto,errflag]=newton(1,10^(-6),40,f,df)
%
%---------------------------------

%---------------------------------------
% PRIMA ITERAZIONE DEL METODO DI NEWTON.
%---------------------------------------
k=1;
x(1,1)=x0;

% VALUTAZIONE FUNZIONE.
fx=feval(funct,x(k,1));

if fx == 0 % x(1) E' SOLUZIONE.
    scarto(k,1)=0;
    errflag=0;
    return;
end

% VALUTAZIONE DERIVATA.
dfx=feval(dfunct,x(k,1));

if dfx == 0 % NEWTON NON PUO' PROSEGUIRE.
    scarto(k)=realmax;
    errflag=1;
    return;
else
    scarto(k,1)=-fx/dfx;
end

fprintf('\n \t [ITER.]: %3.0f',k);
fprintf('[VALORE]: %5.5f [ABS.SCARTO]: %2.2e',x(k,1),abs(scarto(k,1)));

%--------------------------------------------
% ITERAZIONI SUCCESSIVE DEL METODO DI NEWTON.
%--------------------------------------------
while (abs(scarto(k)) > tol) & (k < kmax)
    k=k+1;
    x(k,1)=x(k-1,1)+scarto(k-1,1);

    % VALUTAZIONE FUNZIONE.
    fx=feval(funct,x(k));

    if fx == 0 % x(1) E' SOLUZIONE.
        scarto(k,1)=0;
        errflag=0;
        return;
    end

    % VALUTAZIONE DERIVATA.
    dfx=feval(dfunct,x(k,1));

    if dfx == 0 % NEWTON NON PUO' PROSEGUIRE.
        scarto(k,1)=realmax;
        errflag=1;
        return;
    else
        scarto(k,1)=-fx/dfx;
    end

    fprintf('\n \t [ITER.]: %3.0f',k);
    fprintf('[VALORE]: %5.5f [ABS.SCARTO]: %2.2e',x(k),abs(scarto(k)));
end

% ANALISI DELL'USCITA DI NEWTON.
if (abs(scarto(k,1)) > tol) > 0
    errflag=1;
else
    errflag=0;
end

\end{Verbatim}
\noindent Quindi sul workspace di Matlab/Octave scriviamo
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> [x,k,scarto]=newton(1,10^(-6),20,'f','f1');

 	 [ITER.]:   1[VALORE]: 1.00000 [ABS.SCARTO]: 3.45e-01
 	 [ITER.]:   2[VALORE]: 0.65515 [ABS.SCARTO]: 2.22e-01
 	 [ITER.]:   3[VALORE]: 0.43359 [ABS.SCARTO]: 1.45e-01
 	 [ITER.]:   4[VALORE]: 0.28815 [ABS.SCARTO]: 9.63e-02
 	 [ITER.]:   5[VALORE]: 0.19183 [ABS.SCARTO]: 6.40e-02
 	 [ITER.]:   6[VALORE]: 0.12781 [ABS.SCARTO]: 4.26e-02
 	 [ITER.]:   7[VALORE]: 0.08518 [ABS.SCARTO]: 2.84e-02
 	 [ITER.]:   8[VALORE]: 0.05678 [ABS.SCARTO]: 1.89e-02
 	 [ITER.]:   9[VALORE]: 0.03785 [ABS.SCARTO]: 1.26e-02
 	 [ITER.]:  10[VALORE]: 0.02523 [ABS.SCARTO]: 8.41e-03
 	 [ITER.]:  11[VALORE]: 0.01682 [ABS.SCARTO]: 5.61e-03
 	 [ITER.]:  12[VALORE]: 0.01122 [ABS.SCARTO]: 3.74e-03
 	 [ITER.]:  13[VALORE]: 0.00748 [ABS.SCARTO]: 2.49e-03
 	 [ITER.]:  14[VALORE]: 0.00498 [ABS.SCARTO]: 1.66e-03
 	 [ITER.]:  15[VALORE]: 0.00332 [ABS.SCARTO]: 1.11e-03
 	 [ITER.]:  16[VALORE]: 0.00222 [ABS.SCARTO]: 7.38e-04
 	 [ITER.]:  17[VALORE]: 0.00148 [ABS.SCARTO]: 4.92e-04
 	 [ITER.]:  18[VALORE]: 0.00098 [ABS.SCARTO]: 3.28e-04
 	 [ITER.]:  19[VALORE]: 0.00066 [ABS.SCARTO]: 2.19e-04
 	 [ITER.]:  20[VALORE]: 0.00044 [ABS.SCARTO]: 1.46e-04
>>

\end{Verbatim}
che mostra una convergenza lenta, visto che la radice $x=0$ \mye multipla
(si osservi che $f^{(1)}(0)=\cos{(0)}-1=0$)
\smsp
\begin{remark}
Si osservi che si poteva fare a meno di scrivere le funzioni su file, lavorando direttamente da workspace di Matlab/Octave utilizzando il comando {\tt{inline}}
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> f=inline('sin(x)-x');
>> f1=inline('cos(x)-1');
>> [x,k,scarto]=newton(1,10^(-6),20,'f','f1');

 	 [ITER.]:   1[VALORE]: 1.00000 [ABS.SCARTO]: 3.45e-01
 	 [ITER.]:   2[VALORE]: 0.65515 [ABS.SCARTO]: 2.22e-01
 	 [ITER.]:   3[VALORE]: 0.43359 [ABS.SCARTO]: 1.45e-01
 	 [ITER.]:   4[VALORE]: 0.28815 [ABS.SCARTO]: 9.63e-02
 	 [ITER.]:   5[VALORE]: 0.19183 [ABS.SCARTO]: 6.40e-02
 	 [ITER.]:   6[VALORE]: 0.12781 [ABS.SCARTO]: 4.26e-02
 	 [ITER.]:   7[VALORE]: 0.08518 [ABS.SCARTO]: 2.84e-02
 	 [ITER.]:   8[VALORE]: 0.05678 [ABS.SCARTO]: 1.89e-02
 	 [ITER.]:   9[VALORE]: 0.03785 [ABS.SCARTO]: 1.26e-02
 	 [ITER.]:  10[VALORE]: 0.02523 [ABS.SCARTO]: 8.41e-03
 	 [ITER.]:  11[VALORE]: 0.01682 [ABS.SCARTO]: 5.61e-03
 	 [ITER.]:  12[VALORE]: 0.01122 [ABS.SCARTO]: 3.74e-03
 	 [ITER.]:  13[VALORE]: 0.00748 [ABS.SCARTO]: 2.49e-03
 	 [ITER.]:  14[VALORE]: 0.00498 [ABS.SCARTO]: 1.66e-03
 	 [ITER.]:  15[VALORE]: 0.00332 [ABS.SCARTO]: 1.11e-03
 	 [ITER.]:  16[VALORE]: 0.00222 [ABS.SCARTO]: 7.38e-04
 	 [ITER.]:  17[VALORE]: 0.00148 [ABS.SCARTO]: 4.92e-04
 	 [ITER.]:  18[VALORE]: 0.00098 [ABS.SCARTO]: 3.28e-04
 	 [ITER.]:  19[VALORE]: 0.00066 [ABS.SCARTO]: 2.19e-04
 	 [ITER.]:  20[VALORE]: 0.00044 [ABS.SCARTO]: 1.46e-04
	
>>

\end{Verbatim}
\end{remark}

\subsection{Esercizio sul calcolo delle radici quadrate}

Il calcolo della radice quadrata di un numero reale $a \geq 0$ pu\'o
essere visto come il calcolo dell'unico zero positivo dell'equazione
\begin{equation}\label{rq}
x^{2}-a=0.
\end{equation}
Posto $f(x)=x^{2}-a$, essendo $f^{(1)}(x)=2x$ da (\ref{newton})
otteniamo
\begin{equation}\label{met1}
x_{k+1}=x_{k}+h_{k}, \, \, h_{k}=-
{\frac{x_{k}^{2}-a}{2 \, x_{k}}}
=-{\frac{x_{k}}{2}}+{\frac{a}{2 \, x_{k}}}, \,\, k=0,1,\ldots
\end{equation}
Scriviamo in files {\tt{rad2.m}}, {\tt{drad2.m}}, {\tt{newton.m}} le
seguenti funzioni:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function [fx]=rad2(x)
fx=x.^2-2;

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

function [fx]=drad2(x)
fx=2*x;

\end{Verbatim}
\noindent
A questo punto scriviamo sul workspace di Matlab/Octave
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

[x,k,scarto]=newton(1.41,1e-8,50,'rad2','drad2');

\end{Verbatim}

\noindent
Il metodo stampa sulla stessa workspace
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

 [ITER.]:   1 [VALORE]: 1.41000 [ABS.SCARTO]: 4.22e-003
 [ITER.]:   2 [VALORE]: 1.41422 [ABS.SCARTO]: 6.30e-006
 [ITER.]:   3 [VALORE]: 1.41421 [ABS.SCARTO]: 1.40e-011

\end{Verbatim}
\noindent
e quindi fornisce un'approssimazione della soluzione in sole $3$
iterazioni se partiamo dal punto iniziale $x_0=1.41$.

{\vspace{0.3cm}}

\noindent
{\bf{Nota}}. Osserviamo che
\smsp
\begin{enumerate}
\item posto $x^{\ast}:=\sqrt 2$; \smsp
\item scelti $a$ e $b$ cosicch\mye $0<a<x^{\ast}:=\sqrt 2$, $b > x^{\ast}$; \smsp
\item $f^{\prime}(x)=2\,x \geq 0$, $f^{(2)}(x)=2$ per $x\in [a,b]$. \smsp
\end{enumerate}

Di conseguenza \mye possibile applicare il teorema {\ref{convatk}}, e si conseguenza se $x_0=b$ il metodo di Newton converge a $x^{\ast}=\sqrt(2)$.

In altre parole, scelto arbitrariamente $x_0 > \sqrt 2$ il metodo di Newton converge alla soluzione $\sqrt 2$.


\noindent
Da qui si vede con un po' di tecnica che la successione $x_{k}$ converge a $x_{\ast}$.

\begin{remark}
Si osservi che il metodo di Newton permette di approssimare la radice quadrata di un numero nonnegativo, usando solamente le quattro operazioni $+$, $-$, $*$ e $\backslash$.
\end{remark}

\subsection{Confronto tra il metodo di bisezione e di Newton}
Finora abbiamo implementato mediante una function il metodo di bisezione, usando come
criteri di arresto la semilunghezza dell'intervallo, l'approssimazione
mediante rapporto incrementale del residuo pesato e il numero massimo di
iterazioni.

Abbiamo successivamente eseguito degli esperimenti per esaminare gli zeri
delle funzioni
\smsp
\begin{enumerate}
\item $f(x)=x^2-2$ con $x_1 \in [1,2] $; \smsp
\item $f(x)=\sin (x)-x=0$ \smsp
\end{enumerate}
con particolari tolleranze richieste e numero massimo
di iterazioni.


 \begin{figure}\label{tangenti}
 \centering
   {\includegraphics[scale=0.5,clip]{tangenti.eps}}
 \caption{Grafico che illustra geometricamente le iterazioni del metodo Newton per il calcolo della radice quadrata di 2.}
 \end{figure}


Plottiamo ora, in un unico grafico semi-logaritmico, l'andamento
dell'errore relativo, e su workspace determiniamo la grandezza del residuo pesato relativo e della semi-lunghezza
dell'intervallo relativa per il metodo di bisezione e dell'errore relativo
e dello scarto relativo per il metodo di Newton applicati alla funzione
$$f(x)=x^2-2$$
(con punto iniziale $x_1=1$ per il metodo di Newton), con un numero di
iterazioni pari a $k=50$.

\subsubsection{Commenti all'esercitazione}

 \begin{figure}\label{fig1}
 \centering
   {\includegraphics[scale=0.5,clip]{conv_bn.eps}}
 \caption{Grafico che illustra l'errore relativo (rispetto la soluzione esatta $\sqrt{2}$) del metodo di bisezione e
Newton, rappresentate rispettivamente da {\tt{o}} e {\tt{+}}.}
 \end{figure}

\smsp

\begin{enumerate}
\item Per far eseguire ai programmi esattamente $50$ iterazioni, si pu\myo
impostare il numero massimo di iterazioni a $50$ e la tolleranza a $0$.
\smsp
\item Per disegnare il grafico degli errori, si esegua prima il codice
relativo al metodo di Newton, avendo cura di usare nomi diversi per le due
successioni generate. Conseguentemente, se {\tt{x}} \mye il vettore
contenente la successione prodotta dal metodo di bisezione e {\tt{y}}
quella prodotta dal metodo di Newton
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> [x,k,semilunghezza,residuopesato]=bisezione(1,2,0,0,50,'rad2');

 [IT]: 1 [c]:1.50000 [z1]:1.000 [AMP]:2.50e-001 [WRES]:1.00e-001
 [IT]: 2 [c]:1.25000 [z1]:1.250 [AMP]:1.25e-001 [WRES]:1.59e-001
 [IT]: 3 [c]:1.37500 [z1]:1.375 [AMP]:6.25e-002 [WRES]:3.80e-002
 [IT]: 4 [c]:1.43750 [z1]:1.375 [AMP]:3.13e-002 [WRES]:2.36e-002
 [IT]: 5 [c]:1.40625 [z1]:1.406 [AMP]:1.56e-002 [WRES]:7.90e-003
 [IT]: 6 [c]:1.42188 [z1]:1.406 [AMP]:7.81e-003 [WRES]:7.68e-003
  ....
 [IT]:49 [c]:1.41421 [z1]:1.414 [AMP]:8.88e-016 [WRES]:1.00e-015
 [IT]:50 [c]:1.41421 [z1]:1.414 [AMP]:4.44e-016 [WRES]:1.48e-016
 	
>> [y,kn,scarto]=newton(1,0,50,'rad2','drad2');

 [ITER]:   1 [VALUE]: 1.00000 [STEP]: 5.00e-001
 [ITER]:   2 [VALUE]: 1.50000 [STEP]: 8.33e-002
 [ITER]:   3 [VALUE]: 1.41667 [STEP]: 2.45e-003
 [ITER]:   4 [VALUE]: 1.41422 [STEP]: 2.12e-006
 [ITER]:   5 [VALUE]: 1.41421 [STEP]: 1.59e-012
 [ITER]:   6 [VALUE]: 1.41421 [STEP]: 1.57e-016
  ....
 [ITER]:  50 [VALUE]: 1.41421 [STEP]: 1.57e-016
 	
>> semilogy(abs(x-sqrt(2))/sqrt(2),'r-o'); hold on;
>> semilogy(abs(y-sqrt(2))/sqrt(2),'k-+'); hold on;

\end{Verbatim}
\noindent
produce il grafico degli errori relativi dei due metodi. In rosso a
tondini viene rappresentato l'errore assoluto della
bisezione, in nero coi $+$ quello di Newton.
\end{enumerate}



\section{Esercizi per casa} Calcolare con i metodi di bisezione e Newton un'approssimazione dello zero $\alpha$ di
\begin{eqnarray}
&& 5x = {\exp{(-x)}}  \\
&& x^4-2 = 0 \\
&& x \log{(x)} -1 =0\\
&& (x-1) \log{(x)} =0 \\
&& m  = x - E {\sin{(x)}}, \,\, m=0.8, \, E=2
\end{eqnarray}
Verificare che per le approssimazioni $x^{\ast} \approx\alpha$ fornite sia effettivamente $f({x^{\ast}}) \approx 0$. Posto come $x^{\ast}$ il risultato finale, si disegni il grafico degli errori $|x^{(k)}-x^{\ast}|$ in scala semilogaritmica. Da questa dedurre, se possibile, l'ordine di convergenza del metodo e giustificare teoricamente la risposta.
\section{Facoltativo: Metodo di punto fisso}
Si supponga di dover calcolare $\alpha$ soluzione del problema $x=g(x)$ con $x \in [a,b]$ e $\phi$ sufficientemente derivabile in $[a,b]$. Si dimostra che
\smsp
\begin{theorem} {\label{BAN}}
Si assuma che $g$ sia derivabile con continuit\mya in [a,b], che $g([a,b]) \subseteq [a,b]$ e che se $x \in [a,b]$ sia
$$|g^{\prime}(x)| \leq \lambda < 1.$$
Allora (cf. {\cite[p.80]{ATK}})
\begin{enumerate}
\item $x=g(x)$ ha un'unica soluzione $\alpha$ in $[a,b]$;
\item per qualsiasi scelta di $x_0 \in [a,b]$, posto $$x_{n+1}=g(x_n)$$ \noindent la successione converge ad $\alpha$;
\item si ha $$|\alpha - x_n| \leq {\frac{\lambda^n}{1-\lambda}}|x_1-x_0|.$$
\end{enumerate}
\end{theorem}
\smsp
Implementiamo la successione di punto fisso utilizzando il criterio di arresto dello {\sl{step}} e approssimiamo con tale metodo la soluzione di $$x={\frac{\exp{(-x)}}{5}}$$ con una tolleranza di $10^{(-6)}$ partendo da $x_0 \in [0,1]$ (a piacere).
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function [x,k,scarto]=puntofisso(x0,tol,kmax,g)

% IL METODO CERCA DI APPROSSIMARE LA SOLUZIONE
% DEL PROBLEMA DI PUNTO FISSO "x=g(x)".

% PRIMA ITERAZIONE DEL METODO DI PUNTO FISSO x=g(x).
k=1;
x(1)=x0;
gx=feval(g,x(k));
scarto(k)=abs(gx-x(1));
fprintf('\n \t [ITER.]: %3.0f',k);
fprintf('[VALORE]: %5.5f [ABS.SCARTO]: %2.2e',x(k),abs(scarto(k)));

% ITERAZIONI SUCCESSIVE DEL METODO DI PUNTO FISSO.
while (abs(scarto(k)) > tol) & (k < kmax)
    k=k+1;
    x(k)=gx;
    gx=feval(g,x(k));
    scarto(k)=abs(gx-x(k));
    fprintf('\n \t [ITER.]: %3.0f',k);
    fprintf('[VALORE]: %5.5f [ABS.SCARTO]: %2.2e',x(k),abs(scarto(k)));
end

fprintf('\n');


\end{Verbatim}

\smsp
Testiamo il nostro codice sull'esempio $x=\exp{(-x)}/5$

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

>> % DEFINISCO LA FUNZIONE COME INLINE.
>> h=inline('exp(-x)/5');
>>
>> % CALCOLO L'APPROSSIMAZIONE DELLA SOLUZIONE DI x=h(x).
>> [x,k,scarto]=puntofisso(0,10^(-6),100,h);

 	 [ITER.]:   1[VALORE]: 0.00000 [ABS.SCARTO]: 2.00e-01
 	 [ITER.]:   2[VALORE]: 0.20000 [ABS.SCARTO]: 3.63e-02
 	 [ITER.]:   3[VALORE]: 0.16375 [ABS.SCARTO]: 6.05e-03
 	 [ITER.]:   4[VALORE]: 0.16979 [ABS.SCARTO]: 1.02e-03
 	 [ITER.]:   5[VALORE]: 0.16877 [ABS.SCARTO]: 1.73e-04
 	 [ITER.]:   6[VALORE]: 0.16894 [ABS.SCARTO]: 2.92e-05
 	 [ITER.]:   7[VALORE]: 0.16891 [ABS.SCARTO]: 4.93e-06
 	 [ITER.]:   8[VALORE]: 0.16892 [ABS.SCARTO]: 8.33e-07
>>
>> % VERIFICO CHE x=h(x)
>> format long
>> x(k)

ans =

   0.168916686007286

>> feval(h,x(k))

ans =

   0.168915853145140

>> abs(feval(h,x(k))-x(k))

ans =

     8.328621461384245e-07

>>

\end{Verbatim}

\subsubsection{Esercizio}
Cosa si pu\myo dire relativamente alla sua convergenza? Per mostrare la convergenza del metodo di punto fisso, si pu\myo utilizzare il teorema {\ref{BAN}} per $x_0 \in [0,1]$?

\subsubsection{Esercizio svolto}
Risolvere col metodo di punto fisso l'equazione di Colebrook
$$\sqrt{x}={\frac{1}{-2 \log_{10}{( {\frac{(\epsilon/D)}{3.7}}  + {\frac{2.51}    {{\mbox{Re}}_D  }     }      \cdot {\frac{1}{\sqrt{x}}}           )}}}$$
\noindent in cui il numero di Reynolds \mye ${\mbox{Re}}_D=10^{5}$ ed inoltre $\epsilon/D=0.02$. La soluzione \mye $x \approx 0.0490$.

\noindent
Per prima cosa risolviamo l'equazione $t=C(t)$ dove
$$C(t)={\frac{1}{-2 \log_{10}{( {\frac{(\epsilon/D)}{3.7}}  + {\frac{2.51}    {{\mbox{Re}}_D }     }    \cdot {\frac{1}{t}}           )}}}$$
e una volta ottenuto $t^{\ast}$ tale che $t^{\ast}=C(t^{\ast})$ poniamo $x^{\ast}=(t^{\ast})^2$.

\noindent
Dopo aver effettuato
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> C=inline('-1./(2*log10(0.02/3.7+2.51./(10^5*x)))')

C =

     Inline function:
     C(x) = -1./(2*log10(0.02/3.7+2.51./(10^5*x)))

>> tt=0:0.01:10;
>> yy=feval(C,tt);
>> % GRAFICO DI "t-C(t)"
>> zz=tt-yy;
>> plot(tt,zz)

\end{Verbatim}
deduciamo che la soluzione $t^{\ast}$ sta nell'intervallo ($0,1$) come si vede facilmente dal grafico.

\begin{figure}\label{fig1cole}
\centering
 {\includegraphics[scale=0.5,clip]{colebrook.eps}}
 \caption{Grafico della funzione $t-C(t)$ che mostra che $t-C(t)=0$ per $t\approx 0.2$.}
\end{figure}

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

>> format long
>> [t,k,scarto]=puntofisso(0.2,10^(-13),100,C);

 	 [ITER.]:   1[VALORE]: 0.20000 [ABS.SCARTO]: 2.15e-02
 	 [ITER.]:   2[VALORE]: 0.22151 [ABS.SCARTO]: 9.40e-05
 	 [ITER.]:   3[VALORE]: 0.22142 [ABS.SCARTO]: 3.71e-07
 	 [ITER.]:   4[VALORE]: 0.22142 [ABS.SCARTO]: 1.47e-09
 	 [ITER.]:   5[VALORE]: 0.22142 [ABS.SCARTO]: 5.79e-12
 	 [ITER.]:   6[VALORE]: 0.22142 [ABS.SCARTO]: 2.29e-14

>> xsol=(t(k))^2

xsol =

   0.049026548280573

>>

\end{Verbatim}

\section{Facoltativo: Metodo delle secanti}
Uno dei metodi pi\myu comunemente utilizzati per la risoluzione di equazioni nonlineari \mye quello delle secanti, che nel caso di sistemi di equazioni nonlineari porta (non banalmente) al molto noto metodo di Broyden.

Il metodo delle secanti \mye definito dalla successione

\begin{equation}
x_{n+1}=x_{n}-f(x_{n}) \cdot {  \frac {x_n-x_{n-1}} {f(x_n)-f(x_{n-1})}  }
\end{equation}

\noindent
ottenuta dal metodo di Newton sostituendo $f^{\prime}(x_n)$ col rapporto incrementale
$${  \frac {f(x_n)-f(x_{n-1})} {x_n-x_{n-1}} }.$$

\noindent
Si nota subito che a differenza del metodo di Newton, quello delle secanti richiede due punti iniziali $x_0$, $x_1$ e non necessita ovviamente del calcolo della derivata $f^{\prime}$.

Si dimostra che
\smsp
\begin{theorem}
Se $f \in C^2([a,b])$ con la radice $\alpha \in [a,b]$ e $f^{\prime}(\alpha)$ allora se $x_0$, $x_1$ sono sufficientemente vicine ad $\alpha$ le iterate del metodo delle secanti convergono ad $\alpha$, con ordine di convergenza $$p=  \frac {(1+\sqrt{5})} {2} \approx 1.62.$$
\end{theorem}
\smsp
\noindent
{\bf{Esercizio}}. Aiutandosi con quanto descritto nella routine {\tt{newton.m}}, si implementi il metodo delle secanti e lo si usi per calcolare la radice quadrata di $2$ con una tolleranza di $10^{(-15)}$. Quale criterio di arresto si usi quello step $$|x_{n+1}-x_{n}| \leq {\mbox{toll}}.$$
Se eseguito correttamente si ottiene qualcosa del tipo:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> [x,k,scarto]=secants(1,2,10^(-15),100,'rad2');

 	 [ITER.]:   1 [VALORE]: 2.00000 [ABS.SCARTO]: 1.00e+000
 	 [ITER.]:   2 [VALORE]: 1.33333 [ABS.SCARTO]: 6.67e-001
 	 [ITER.]:   3 [VALORE]: 1.40000 [ABS.SCARTO]: 6.67e-002
 	 [ITER.]:   4 [VALORE]: 1.41463 [ABS.SCARTO]: 1.46e-002
 	 [ITER.]:   5 [VALORE]: 1.41421 [ABS.SCARTO]: 4.23e-004
 	 [ITER.]:   6 [VALORE]: 1.41421 [ABS.SCARTO]: 2.12e-006
 	 [ITER.]:   7 [VALORE]: 1.41421 [ABS.SCARTO]: 3.16e-010
 	 [ITER.]:   8 [VALORE]: 1.41421 [ABS.SCARTO]: 3.14e-016
 	
\end{Verbatim}
\noindent
supposto che il valore alla prima iterazione sia quello di $x(2)$ che nell'esempio \mye $2$. Calcoliamo una stima della velocit\mya di convergenza
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> abs_scarto=abs(scarto);
>> L=length(abs_scarto);
>> ratios=abs_scarto(2:L)./abs_scarto(1:L-1);
>> L1=length(ratios);
>> p_approx=log(ratios(2:L1))./log(ratios(1:L1-1));
>> format long
>> p_approx
p_approx =
   5.67887358726758
   0.65854134728041
   2.33747974959565
   1.49349105655602
   1.66495845984246
   1.56815727154018

\end{Verbatim}
\noindent
e quindi la velocit\mya di convergenza stimata \mye prossima a quella teorica $(1+\sqrt{5})/2 \approx 1.62$.

Di seguito si esegua il calcolo della radice cubica del proprio numero di matricola.

\section{Facoltativo: Un'equazione polinomiale}
Si supponga di dover calcolare la pi\myu grande radice dell'equazione
\begin{equation}
x^6-x+1=0
\end{equation}
A tal proposito, aiutandosi con pico o l'editor preferito, poniamo
{\tt{eqtype=1}} in {\tt{g.m}} e {\tt{dg.m}}.


Eseguiamo un plot per avere un'idea della disposizione degli zeri.
Da una regola dovuta a Cartesio, se il polinomio $p$ di cui calcolare gli zeri $z_i$ \mye
\begin{equation}
p(x):=\sum_{j=0}^{n} a_j x^j
\end{equation}
allora
\begin{equation}
|z_i| \leq 1+\max_{0 \leq i \leq n-1} \frac{|a_i|}{|a_n|}
\end{equation}
Nel nostro caso, essendo
$$a_0=1, \, a_1=-1, \, a_6=1, \, a_2=\ldots=a_5=0$$
\noindent
si ha quindi  $|z_i| \leq 2$.

\noindent
Eseguiamo dal workspace di Matlab/Octave il comando
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> x=-2:0.001:2; y=g(x); plot(x,y,'r-'); hold on;
>> z=zeros(size(y)); plot(x,z,'g-'); hold on;

\end{Verbatim}
\noindent Come risultato abbiamo il grafico in figura. Da questo si vede che la radice (positiva) cercata risulta essere un numero nell'intervallo
[$1.0$,$1.5$]. Zoomando pi\myu volte si vede che appartiene all'intervallo [$1.11$,$1.16$]. Lanciamo
quindi il metodo di bisezione (se ne osservi l'applicabilit\mya!) ottenendo

 \begin{figure}\label{fig1b}
 \centering
   {\includegraphics[scale=0.5,clip]{roots_poly.eps}}
   {\includegraphics[scale=0.5,clip]{roots_poly_zoom.eps}}
 \caption{Grafico del polinomio $x^6-x-1$.}
 \end{figure}

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

>> [c,k,semilunghezza,residuopesato]=bisezione(1.12,1.16,10^(-8),10^(-8),50,'g');

 [IT]:   1 [c]: 1.14000 [z1]: 1.120 [AMP]: 1.00e-002 [WRES]: 5.47e-003
 [IT]:   2 [c]: 1.13000 [z1]: 1.130 [AMP]: 5.00e-003 [WRES]: 4.66e-003
 [IT]:   3 [c]: 1.13500 [z1]: 1.130 [AMP]: 2.50e-003 [WRES]: 2.79e-004
 [IT]:   4 [c]: 1.13250 [z1]: 1.133 [AMP]: 1.25e-003 [WRES]: 2.22e-003
 [IT]:   5 [c]: 1.13375 [z1]: 1.134 [AMP]: 6.25e-004 [WRES]: 9.73e-004
 [IT]:   6 [c]: 1.13438 [z1]: 1.134 [AMP]: 3.12e-004 [WRES]: 3.49e-004
 [IT]:   7 [c]: 1.13469 [z1]: 1.135 [AMP]: 1.56e-004 [WRES]: 3.66e-005
 [IT]:   8 [c]: 1.13484 [z1]: 1.135 [AMP]: 7.81e-005 [WRES]: 1.20e-004
 [IT]:   9 [c]: 1.13477 [z1]: 1.135 [AMP]: 3.91e-005 [WRES]: 4.15e-005
 [IT]:  10 [c]: 1.13473 [z1]: 1.135 [AMP]: 1.95e-005 [WRES]: 2.42e-006
 [IT]:  11 [c]: 1.13471 [z1]: 1.135 [AMP]: 9.77e-006 [WRES]: 1.71e-005
 [IT]:  12 [c]: 1.13472 [z1]: 1.135 [AMP]: 4.88e-006 [WRES]: 7.34e-006
 [IT]:  13 [c]: 1.13472 [z1]: 1.135 [AMP]: 2.44e-006 [WRES]: 2.46e-006
 [IT]:  14 [c]: 1.13472 [z1]: 1.135 [AMP]: 1.22e-006 [WRES]: 1.73e-008
 [IT]:  15 [c]: 1.13473 [z1]: 1.135 [AMP]: 6.10e-007 [WRES]: 1.20e-006
 [IT]:  16 [c]: 1.13472 [z1]: 1.135 [AMP]: 3.05e-007 [WRES]: 5.93e-007
 [IT]:  17 [c]: 1.13472 [z1]: 1.135 [AMP]: 1.53e-007 [WRES]: 2.88e-007
 [IT]:  18 [c]: 1.13472 [z1]: 1.135 [AMP]: 7.63e-008 [WRES]: 1.35e-007
 [IT]:  19 [c]: 1.13472 [z1]: 1.135 [AMP]: 3.81e-008 [WRES]: 5.90e-008
 [IT]:  20 [c]: 1.13472 [z1]: 1.135 [AMP]: 1.91e-008 [WRES]: 2.08e-008
 	
>>

\end{Verbatim}
\noindent
Desiderando avere almeno 8 cifre decimali come richiesto dalla tolleranza, dopo aver ricordato che
$c$ \mye un vettore con la {\sl{storia}} dell'approssimazione dello zero cercato, come si evince da
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> c

c =
  Columns 1 through 7
    1.1400    1.1300    1.1350    1.1325    1.1338    1.1344    1.1347
  Columns 8 through 14
    1.1348    1.1348    1.1347    1.1347    1.1347    1.1347    1.1347
  Columns 15 through 21
    1.1347    1.1347    1.1347    1.1347    1.1347    1.1347    1.1347

>>

\end{Verbatim}

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

>> format long; M=length(c); risultato=c(M)

risultato =
   1.13472414016724

>>

\end{Verbatim}

\noindent
Se intendiamo usare il metodo di Newton, bisogna calcolare la derivata di $g(x):=x^6-x-1$ che \mye uguale a $Dg(x):=6\,x^5-1$. Questa \mye implementata
nel file {\tt{dg.m}} e basta settare il parametro {\tt{eqtype=1;}} per poterla valutare.

\noindent Dopo aver fatto questa modifica, testiamo il metodo di Newton con punto iniziale $x_0=1.12$, tolleranza $10^{-8}$ e numero massimo di
iterazioni uguale a $50$, ottenendo:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> [x,k,scarto,errflag]=newton(1.12,1e-8,50,'g','dg');

 [ITER]:   1 [VALUE]: 1.12000 [STEP]: 1.53e-002
 [ITER]:   2 [VALUE]: 1.13527 [STEP]: 5.43e-004
 [ITER]:   3 [VALUE]: 1.13472 [STEP]: 7.14e-007
 [ITER]:   4 [VALUE]: 1.13472 [STEP]: 1.23e-012

>> format long;
>> x

x =
   1.12000000000000   1.13526807498409   1.13472485264634   1.13472413840275

>> M=length(x); risultato=x(M)

risultato =
   1.13472413840275

>>

\end{Verbatim}
\subsection{Facoltativo: Un esempio sul calcolo dello zero di una funzione}
Consideriamo la funzione
$$f_2(x):=x-6.28-\sin(x)$$
e supponiamo di doverne calcolare uno zero $x^{\ast}$. Notiamo che la soluzione non \mye $2\pi$. Infatti
$$f_2(2\pi):=2\pi-6.28-\sin{(2\pi)}\approx 0.00318530717959$$

\noindent
come si vede dal workspace di Matlab/Octave
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> format long
>> 2*pi-6.28-sin(2*pi)
ans =
   0.00318530717959
>>

\end{Verbatim}
Essendo $\sin(x) \in [-1,1]$ necessariamente $6.28+\sin(x)$ sta nell'intervallo [$5.28,7.28$] e
visto che
$$x-6.28-\sin(x)=0 \Leftrightarrow x=6.28+\sin(x)$$
deduciamo che $x^{\ast} \in [5.28,7.28]$. Eseguiamo quindi un plot in tale intervallo. Dopo aver settato {\tt{eqtype=2}} tanto in {\tt{g}} quanto in
{\tt{dg}}, digitiamo nel workspace di Matlab/Octave
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> x=5.28:0.001:7.28; y=g(x); plot(x,y); hold on;
>> x=5.28:0.001:7.28; y=zeros(size(x)); plot(x,y,'g-'); hold on;

\end{Verbatim}
\noindent (si confronti con la ralativa figura). Dopo qualche zoom intorno allo zero della funzione, si evince che lo zero cercato \mye nell'intervallo
[$6,6.04$].

 \begin{figure}\label{fig1c}
 \centering
   {\includegraphics[scale=0.5,clip]{roots_f2_wide.eps}}
   {\includegraphics[scale=0.5,clip]{roots_f2.eps}}
   {\includegraphics[scale=0.5,clip]{roots_f2_zoom.eps}}
 \caption{Grafico del polinomio $x-6.28-\sin(x)$ su varie scale, vicino alla radice $x^{\ast} \approx 6.01$.}
 \end{figure}

\noindent Dal grafico nell'intervallo pi\myu ampio, si vede che il metodo delle tangenti (vederlo geometricamente!) converge partendo da $x_0=4$.
Otteniamo quale risultato
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> [x,k,scarto,errflag]=newton(4,1e-8,50,'g','dg');

 [ITER]:   1 [VALUE]: 4.00000 [STEP]: 9.21e-001
 [ITER]:   2 [VALUE]: 4.92112 [STEP]: 4.80e-001
 [ITER]:   3 [VALUE]: 5.40118 [STEP]: 2.93e-001
 [ITER]:   4 [VALUE]: 5.69428 [STEP]: 1.80e-001
 [ITER]:   5 [VALUE]: 5.87397 [STEP]: 9.86e-002
 [ITER]:   6 [VALUE]: 5.97256 [STEP]: 3.73e-002
 [ITER]:   7 [VALUE]: 6.00988 [STEP]: 5.51e-003
 [ITER]:   8 [VALUE]: 6.01539 [STEP]: 1.14e-004
 [ITER]:   9 [VALUE]: 6.01550 [STEP]: 4.85e-008
 [ITER]:  10 [VALUE]: 6.01550 [STEP]: 7.79e-015

>>

\end{Verbatim}
\noindent il che mostra la convergenza piuttosto lenta del metodo. Similmente, il metodo converge pure partendo da $x_0=8$, ma nuovamente i risultati non sono entusiasmanti
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> [x,k,scarto,errflag]=newton(8,1e-8,50,'g','dg');

 [ITER]:   1 [VALUE]: 8.00000 [STEP]: 6.38e-001
 [ITER]:   2 [VALUE]: 7.36216 [STEP]: 3.80e-001
 [ITER]:   3 [VALUE]: 6.98191 [STEP]: 2.50e-001
 [ITER]:   4 [VALUE]: 6.73155 [STEP]: 1.83e-001
 [ITER]:   5 [VALUE]: 6.54886 [STEP]: 1.80e-001
 [ITER]:   6 [VALUE]: 6.36930 [STEP]: 8.88e-001
 [ITER]:   7 [VALUE]: 5.48106 [STEP]: 2.63e-001
 [ITER]:   8 [VALUE]: 5.74386 [STEP]: 1.59e-001
 [ITER]:   9 [VALUE]: 5.90295 [STEP]: 8.28e-002
 [ITER]:  10 [VALUE]: 5.98571 [STEP]: 2.69e-002
 [ITER]:  11 [VALUE]: 6.01264 [STEP]: 2.84e-003
 [ITER]:  12 [VALUE]: 6.01547 [STEP]: 3.01e-005
 [ITER]:  13 [VALUE]: 6.01550 [STEP]: 3.35e-009

>>

\end{Verbatim}
\noindent
Dalle ultime iterazioni si capisce che la convergenza \mye comunque superlineare.

\subsection{Facoltativo. Zeri multipli di una funzione e metodo di Newton: un esempio}

Abbiamo visto che il metodo di Newton non ha convergenza
quadratica nel caso del problema $f(x)=0$ in cui $x^{\ast}$ sia uno zero multiplo di molteplicit\mya $p > 1$ cio\mye tale che
$$f(x^{\ast})=f^{(1)}(x^{\ast})=\ldots=f^{(p-1)}(x^{\ast})=0.$$
Per esempio, si calcolino con il metodo di Newton gli zeri multipli della funzione
$$f_3(x):=\exp(x^2)-1$$
implementata in {\tt{g.m}} con {\tt{eqtype=3}} (ricordarsi di settare pure {\tt{dg.m}}!).
La funzione ha uno zero multiplo in $x^{\ast}=0$, essendo
$$Df_3(x):=2x\exp(x^2), \, \, f_3(x^{\ast})=Df_3(x^{\ast})=0.$$
D'altra parte
$$f_3(x):=\exp(x^2)-1=0 \Leftrightarrow \exp(x^2)=1 \Leftrightarrow x^2=\log{({\exp(x^2)})}=\log{({1})}=0$$
il he nuovamente mostra che la radice $x=0$ \mye multipla (calcolare la derivata prima di $f$ nella radice $x=0$).

Dal plot dell'errore del metodo in scala semi-logaritmica, valutare l'ordine di convergenza del metodo.

{\bf{Risoluzione.}}
Settando {\tt{eqtype=3}} in {\tt{g.m}} e {\tt{dg.m}}, lanciamo il metodo di Newton per diversi valori iniziali.
Il metodo non \mye globalmente convergente. Ad esempio:
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> [x,k,scarto,errflag]=newton(8,1e-8,50,'g','dg');

 [ITER]:   1 [VALUE]: 8.00000 [STEP]: 5.44e+027
 [ITER]:   2 [VALUE]: -5443167958782763500000000000.00000 [STEP]: Inf
 [ITER]:   3 [VALUE]:  -Inf [STEP]: NaN

>>

\end{Verbatim}
\noindent Una radice \mye naturalmente $x^{\ast}=0$, ed essendo $\exp(x^2) > 1$ se $x \neq 0$, essa \mye pure l'unica.
Come visto \mye doppia (ha cio\mye molteplicit\mya 2), e dalla teoria si pu\myo
provare che di conseguenza la convergenza di Newton \mye lineare. Vediamolo direttamente partendo da $x_0=1$.
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> [x,k,scarto,errflag]=newton(1,1e-8,50,'g','dg');

 [ITER]:   1 [VALUE]: 1.00000 [STEP]: 3.16e-001
 [ITER]:   2 [VALUE]: 0.68394 [STEP]: 2.73e-001
 [ITER]:   3 [VALUE]: 0.41081 [STEP]: 1.89e-001
 [ITER]:   4 [VALUE]: 0.22180 [STEP]: 1.08e-001
 [ITER]:   5 [VALUE]: 0.11359 [STEP]: 5.64e-002
 [ITER]:   6 [VALUE]: 0.05716 [STEP]: 2.85e-002
 [ITER]:   7 [VALUE]: 0.02863 [STEP]: 1.43e-002
 [ITER]:   8 [VALUE]: 0.01432 [STEP]: 7.16e-003
 [ITER]:   9 [VALUE]: 0.00716 [STEP]: 3.58e-003
 [ITER]:  10 [VALUE]: 0.00358 [STEP]: 1.79e-003
 [ITER]:  11 [VALUE]: 0.00179 [STEP]: 8.95e-004
 [ITER]:  12 [VALUE]: 0.00090 [STEP]: 4.48e-004
 [ITER]:  13 [VALUE]: 0.00045 [STEP]: 2.24e-004
 [ITER]:  14 [VALUE]: 0.00022 [STEP]: 1.12e-004
 [ITER]:  15 [VALUE]: 0.00011 [STEP]: 5.59e-005
 [ITER]:  16 [VALUE]: 0.00006 [STEP]: 2.80e-005
 [ITER]:  17 [VALUE]: 0.00003 [STEP]: 1.40e-005
 [ITER]:  18 [VALUE]: 0.00001 [STEP]: 6.99e-006
 [ITER]:  19 [VALUE]: 0.00001 [STEP]: 3.50e-006
 [ITER]:  20 [VALUE]: 0.00000 [STEP]: 1.75e-006
 [ITER]:  21 [VALUE]: 0.00000 [STEP]: 8.74e-007
 [ITER]:  22 [VALUE]: 0.00000 [STEP]: 4.37e-007
 [ITER]:  23 [VALUE]: 0.00000 [STEP]: 2.18e-007
 [ITER]:  24 [VALUE]: 0.00000 [STEP]: 1.09e-007
 [ITER]:  25 [VALUE]: 0.00000 [STEP]: 5.48e-008
 [ITER]:  26 [VALUE]: 0.00000 [STEP]: 2.64e-008
 [ITER]:  27 [VALUE]: 0.00000 [STEP]: 1.57e-008
 [ITER]:  28 [VALUE]: 0.00000 [STEP]: 8.85e-009
 	
>>

\end{Verbatim}
Vediamo direttamente l'ordine di convergenza. Visto che la soluzione \mye $x^{\ast}=0$ necessariamente
$$ e_k=|x^{(k)}-x^{\ast}|=|x^{(k)}|.   $$
Se il metodo ha ordine $p$, per $k$ sufficientemente grande, esiste C indipendente da $k$ tale che
$$ e_{k+1} \approx C e_{k}^p   $$
$$ e_{k+2} \approx C e_{k+1}^p   $$
e quindi dividendo membro a membro abbiamo
$$ \frac{e_{k+2}}{e_{k+1}}  \approx \left( \frac{e_{k+1}}{e_{k}} \right)^p.$$
Calcolando il logaritmo di ambo i membri
$$ \log \left( \frac{e_{k+2}}{e_{k+1}} \right) \approx \log \left( \frac{e_{k+1}}{e_{k}} \right)^p = p \log \left( \frac{e_{k+1}}{e_{k}} \right)$$
da cui
$$p \approx \frac{\log \left( \frac{e_{k+2}}{e_{k+1}} \right)}{\log \left( \frac{e_{k+1}}{e_{k}} \right)}.$$
Nel nostro caso, essendo $e_{k}=|x^{(k)}|$ abbiamo
$$p \approx \frac{\log \left( \frac{ |x^{(k+2)}| }{ |x^{(k+1)}| } \right)}{\log \left( \frac{ |x^{(k+1)}| }{ |x^{(k)}| } \right)}.$$
Dal workspace di Matlab/Octave
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

>> absx=abs(x);
>> L=length(absx);
>> ratios=absx(2:L)./absx(1:L-1);
>> L1=length(ratios);
>> p_approx=log(ratios(2:L1))./log(ratios(1:L1-1));
>> format long
>> p_approx'

ans =

   1.34180444561870
   1.20915162745740
   1.08581578729537
   1.02616296799749
   1.00694914160602
   1.00176551299514
   1.00044319047015
   1.00011091165740
   1.00002773504230
   1.00000693421959
   1.00000173374970
   1.00000043279671
   1.00000011084991
   1.00000003726416
   0.99999996879110
   1.00000015214145
   0.99999983966276
   1.00000320369014
   0.99999323164192
   1.00004521264807
   0.99976874978756
   0.99962083100483
   0.99886365484392
   1.00329108837037
   0.95044659133176
   1.23127847870172

>>

\end{Verbatim}
da cui si evince sperimentalmente che effettivamente la convergenza \mye lineare.

\subsection{Facoltativo: Function {\tt{fsolve}} }
Si testi la funzione {\tt{fsolve}} di Matlab o GNU Octave sui casi
precedenti, aiutandosi con gli esempi forniti dal comando {\tt{help -i
fsolve}}.

\noindent
In GNU-Octave (con cygwin, su Windows XP) versione $2.1.73$ abbiamo
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

octave:1> [x,info,msg]=fsolve('rad2',1);
warning: in /usr/lib/octave/2.1.73/oct/i686-pc-cygwin/fsolve.oct
near line 14, column 13:
warning: time stamp for `/cygdrive/d/CORSI_MATLAB/MIEI/SM_EQNON
LINEARI_2007/MFILES/rad2.m' is in the future
octave:2> x
x = 1.4142

\end{Verbatim}

\noindent
In Matlab $6.1.0.450$

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

>> fun = inline('x.^2-2');
>> x = fsolve(fun,[1 2],optimset('Display','off'));
>> format long;
>> x
x =
   1.41421356237470   1.41421356237469

\end{Verbatim}
\section{Facoltativo. Nota storica}
Il metodo di Newton \mye anche noto come metodo di Newton-Raphson, in quanto sviluppato da Newton nell'opera {\sl{Analysis Aequationes Numero Terminorum Infinitas}} (1666, scritto nel 1671 e quindi pubblicato nel 1711 da William Jones) per calcolare la radice $2.09455148154233\ldots$ di $x^3-2 \, x-5$. Il metodo forse derivava da alcune idee di Francois Viete che a sua volta si rifaceva al lavoro del matematico persiano Sharaf al-Din al-Tusi. Erone di Alexandria a sua volta utilizzava idee in qualche modo simili.

La questione venne ripresa nel {\sl{De metodis fluxionum et serierum infinitarum}} (scritto nel 1671, tradotto e pubblicato come {\sl{Method of Fluxions}} nel 1736 da John Colson). Il metodo descritto era molto differente
da quello attuale, e veniva sostanzialmente applicato a polinomi {\cite[p. 413]{COM}}, {\cite{WIK4}}.

Il metodo di Newton venne pure pubblicato nel {\sl{A Treatise of Algebra both Historical and Practical}} (1685) di John Wallis.

Nel 1690, Joseph Raphson pubblic\myo nel {\sl{Analysis aequationum universalis}} una versione maggiormente algoritmica, iterativa e meno complicata di quella di Newton. Quindi Simpson, nel 1740, generalizz\myo queste idee per equazioni nonlineari e sistemi nonlineari di piccole dimensioni.

Nel 1768, Jean-Raymond Mourraille nel {\sl{Trait\myer de la r{\myer}solution des {\myer}quation num{\myer}riques}} cominci\myo a discutere la convergenza del metodo di Newton, sottolineando il legame con l'aspetto geometrico delle tangenti.

Pure Fourier diede impulso a questi studi nel {\sl{Question d'analyse algebrique}} (1818).


 \begin{figure}\label{fig1b}
 \centering
 {\includegraphics[scale=0.75,clip]{newton.eps}}
   {\includegraphics[scale=0.7,clip]{fourier.eps}}

 \caption{Ritratti di Newton (1642-1727) e Fourier (1768-1830).}
 \end{figure}

\section{Esercizio facoltativo: il metodo di Halley}
Il metodo di Halley \cite{WIK31} approssima uno zero $x^\ast$ dell'equazione $f(x)=0$ con la sequenza di iterazioni
$$x_{n+1} = x_n - \frac {2 f(x_n) f'(x_n)} {2 {[f'(x_n)]}^2 - f(x_n) f''(x_n)}$$
a partire da un punto iniziale $x_0$. Ovviamente, si richiede l'esistenza della derivata seconda.

Se $f$ \mye differenziabile tre volte con continuit\mya e $x^\ast$ non \mye uno zero delle sue derivate, allora si pu\myo dimostrare (non immediato!) che in un certo intorno di $x^\ast$ le iterate $x_n$ soddisfano

$$| x_{n+1} - x^\ast | \le K \cdot {| x_n - x^\ast |}^3, \text{ per qualche }K > 0.\!$$

\noindent
cio\mye con convergenza cubica.

\noindent
Quale esercizio si
\smsp
\begin{enumerate}
\item approssimi col metodo di Halley lo zero di $$f_2(x):=x-6.28-\sin(x),$$ partendo da $x_0=4$; \smsp
\item approssimi col metodo di Halley la radice quadrata di $2$ partendo da $x_0=1$; \smsp
\item approssimi col metodo di Halley la radice cubica di $2$ partendo da $x_0=1$. \smsp
\end{enumerate}

\noindent
interrompendo il processo quando lo step $|x_{n+1}-x_{n}|$ \mye inferiore di $10^{-14}$.
Si paragoni quindi il metodo di Halley con quello di Newton, sugli stessi problemi: quale dei due sembra convergere pi\myu velocemente?

\section{Alcuni esercizi facoltativi}
\smsp
\begin{enumerate}
\item Scritto un file {\tt{g.m}}
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function y=g(x)

eqtype=1;

switch eqtype
case 1
    y=x.^6-x-1;
case 2
    y=x-6.28-sin(x);
case 3
    y=exp(x.^2)-1;
end

\end{Verbatim}
\noindent
e {\tt{dg.m}} relativo alle corrispondenti derivate
\begin{Verbatim}[fontsize=\small, formatcom=\color{blue}]

function y=dg(x)

eqtype=1;

switch eqtype
case 1
    y=6*(x.^5)-1;
case 2
    y=1-cos(x);
case 3
    y=2*x.*exp(x.^2);
end

\end{Verbatim}
\noindent
si studi l'equazione $g(x)=0$ al variare di {\tt{eqtype}} tra 1 e 3.
\smsp
\item Si risolva numericamente il problema di
calcolare la radice cubica di un numero, col metodo di Newton.
\smsp
\item Si risolva numericamente il problema di
calcolare lo zero di $f(x)={\mbox{arctan}} (x)$, col metodo di Newton partendo prima da $x_0=1.2$ ed una seconda volta da $x_0=1.4$. Convergono in entrambi casi?
\smsp
\item Si calcoli il reciproco di un numero $c$, risolvendo l'equazione $\frac{1}{x} -c =0$ col metodo di Newton. Osservare che questo era utile per vecchi computers che non avevano un chip che faceva la divisione, bens\myi uno che
faceva la moltiplicazione (cf. {\cite[p.35-39]{GAM}}).
\smsp
\item Nota la soluzione del problema $x^2-a=0$
(si ottiene
con il comando
{\tt{sqrt(a)}}) si stimi l'ordine di convergenza del metodo di
Newton (\ref{met1}) per il calcolo della radice quadrata di un numero.
A questo proposito conoscendo $|e_{k+1}|$, $|e_{k}|$, $|e_{k-1}|$ e
ritenendo
costante la costante asintotica di errore $C$ calcolare un'approssimazione
di $p$, via il calcolo di un logaritmo. Si ricordi che
$$p \approx \frac{\log \left( \frac{ |e^{(k+2)}| }{ |e^{(k+1)}| } \right)}{\log \left( \frac{ |e^{(k+1)}| }{ |e^{(k)}| } \right)}.$$

\end{enumerate}


\subsection{Online}

Per approfondimenti, si considerino le pagine web
\smsp
\begin{enumerate}
\item \href{http://it.wikipedia.org/wiki/Calcolo\_dello\_zero\_di\_una\_funzione}{http://it.wikipedia.org/wiki/Calcolo\_dello\_zero\_di\_una\_funzione} \smsp
\item \href{http://it.wikipedia.org/wiki/Metodo\_della\_bisezione}{http://it.wikipedia.org/wiki/Metodo\_della\_bisezione}\smsp
\item \href{http://it.wikipedia.org/wiki/Metodo\_delle\_tangenti}{http://it.wikipedia.org/wiki/Metodo\_delle\_tangenti}\smsp
\item \href{http://en.wikipedia.org/wiki/Newton's\_method}{http://en.wikipedia.org/wiki/Newton's\_method}\smsp
\item \href{http://en.wikipedia.org/wiki/Joseph\_Raphson}{http://en.wikipedia.org/wiki/Joseph\_Raphson}\smsp
\item \href{http://www-history.mcs.st-andrews.ac.uk/Biographies/Fourier.html}{http://www-history.mcs.st-andrews.ac.uk/Biographies/Fourier.html}\smsp
\item \href{http://www-history.mcs.st-andrews.ac.uk/Biographies/Newton.html}{http://www-history.mcs.st-andrews.ac.uk/Biographies/Newton.html}\smsp
\item \href{http://www-history.mcs.st-andrews.ac.uk/Biographies/Raphson.html}{http://www-history.mcs.st-andrews.ac.uk/Biographies/Raphson.html}\smsp
\end{enumerate}

\section{Frasi celebri}
\smsp
\begin{enumerate}
\item Mathematical Analysis is as extensive as nature herself. (Fourier) \smsp
\item Laplace seems quite young; his voice is quiet but clear, and he speaks precisely, though not very fluently; his appearance is pleasant, and he dresses very simply; he is of medium height. His teaching of mathematics is in no way remarkable and he covers the material very rapidly. (Fourier) \smsp
\item Monge has a loud voice and he is energetic, ingenious and very learned. It is well known that his talent is particularly for geometry, physics and chemistry. The subject that he teaches is a fascinating one, and he describes it with the greatest possible clarity. He is even considered to be too clear, or, rather to deal with his material too slowly. He gives individual practical lessons to his students. He speaks colloquially, and for the most part precisely. He is not only to be considered for his great knowledge but is also greatly respected in public and in private. His appearance is very ordinary. (Fourier) \smsp
\item Lagrange, the foremost scholar of Europe, appears to be between 50 and 60 years old, though he is in fact younger; he has a strong Italian accent and pronounces an 's' as if it were a 'z'; he dresses very quietly, in black or brown; he speaks colloquially and with some difficulty, with the hesitant simplicity of a child. Everyone knows that he is an extraordinary man, but one needs to have seen him to recognise him as a great one. He speaks only in discussion, and some of what he says excites ridicule. The other day he said "There are a lot of important things to be said on this subject, but I shall not say them". The students are incapable of appreciating his genius, but the teachers make up for that. (Fourier) \smsp
\item In the absence of any other proof, the thumb alone would convince me of God's existence. (Newton) \smsp
\item Truth is ever to be found in the simplicity, and not in the multiplicity and confusion of things. (Newton) \smsp
\item Tact is the art of making a point without making an enemy. (Newton) \smsp
\item We build too many walls and not enough bridges. (Newton) \smsp
\item Oh Diamond! Diamond! Thou little knowest the mischief done! (Said to a pet dog who knocked over a candle and set fire to his papers). (Newton) \smsp
\item Gravitation is not responsible for people falling in love. (Newton?) \smsp
\end{enumerate}

\begin{thebibliography}{99}
\bibitem{ATK} K. Atkinson, {\em An Introduction to Numerical Analysis\/},
Wiley, (1989).

\bibitem{COM} V. Comincioli, {\em Analisi Numerica, metodi modelli
applicazioni\/}, McGraw-Hill, (1990).


\bibitem{DSC} J.E. Dennis, Jr. e R.B. Schnabel, {\em Numerical Methods for Unconstrained Optimization and Nonlinear Equations\/}, SIAM, (1996).
\bibitem{GAM} G. Gambolati, {\em Elementi di Calcolo Numerico\/}, Edizioni Libreria Cortina, (1988).
\bibitem{KEL} C.T. Kelley, {\em Iterative Methods for Linear and Nonlinear Equations\/}, SIAM Frontiers in Applied Mathematics. SIAM, Philadelphia, (1995).
\bibitem{QUA} A. Quarteroni, F. Saleri {\em Introduzione al Calcolo
Scientifico\/}, Springer, (2002).

\bibitem{MTAF} Mac Tutor (Fourier), \\  {\href{http://www-history.mcs.st-andrews.ac.uk/Biographies/Fourier.html}{http://www-history.mcs.st-andrews.ac.uk/Biographies/Fourier.html}}.
\bibitem{MTAN} Mac Tutor (Newton), \\   {\href{http://www-history.mcs.st-andrews.ac.uk/Biographies/Newton.html}{http://www-history.mcs.st-andrews.ac.uk/Biographies/Newton.html}}.
\bibitem{MTAR} Mac Tutor (Raphson), \\   {\href{http://www-history.mcs.st-andrews.ac.uk/Biographies/Raphson.html}{http://www-history.mcs.st-andrews.ac.uk/Biographies/Raphson.html}}.
\bibitem{ROD} G. Rodriguez, {\em Algoritmi Numerici\/}, Pitagora Editrice, (2008).
\bibitem{Wikipedia} Wikipedia (Newton), \\   {\href{http://en.wikiquote.org/wiki/Newton}{http://en.wikiquote.org/wiki/Newton}}.
\bibitem{WIK1} Wikipedia (Calcolo dello zero di una funzione), \\   {\href{http://it.wikipedia.org/wiki/Calcolo\_dello\_zero\_di\_una\_funzione}{http://it.wikipedia.org/wiki/Calcolo\_dello\_zero\_di\_una\_funzione}}.
\bibitem{WIK2} Wikipedia (Metodo della bisezione), \\   {\href{http://it.wikipedia.org/wiki/Metodo\_della\_bisezione}{http://it.wikipedia.org/wiki/Metodo\_della\_bisezione}}.
\bibitem{WIK3} Wikipedia (Metodo delle tangenti), \\   {\href{http://it.wikipedia.org/wiki/Metodo\_delle\_tangenti}{http://it.wikipedia.org/wiki/Metodo\_delle\_tangenti}}.
\bibitem{WIK31} Wikipedia (Halley's method), \\
{\href{http://en.wikipedia.org/wiki/Halley's\_method}{http://en.wikipedia.org/wiki/Halley's\_method}}.
\bibitem{WIK4} Wikipedia (Newton's method), \\   {\href{http://en.wikipedia.org/wiki/Newton's\_method}{http://en.wikipedia.org/wiki/Newton's\_method}}.
\bibitem{WIK5} Wikipedia (Joseph Raphson), \\   {\href{http://en.wikipedia.org/wiki/Joseph\_Raphson}{http://en.wikipedia.org/wiki/Joseph\_Raphson}}.

\bibitem{ZEI} E. Zeidler, {\em Nonlinear Functional Analysis and its Applications: Part 1: Fixed-Point Theorems\/}, Springer, (1998).

\end{thebibliography}

\end{document}
