% -*- mode: Noweb; noweb-code-mode: c-mode -*-

\hypertarget{sec:newton}{\section{Numerical Analysis}}
\phantomsection
\label{sec:newton:label}

\begin{fullwidth}

\epigraph{An approximate answer to the right problem is worth
a good deal more than an exact answer to an approximate problem.}
{John Tukey}

Much of early computing revolved around numerical algorithms 
that produce values that approximate the \true value.
Approximation algorithms can be terminated once a 
computed value is determined to be \emph{good enough}.
To illustrate numerical algorithms, let's develop Newton's method for
computing an $x$ such that $f(x)=0$.
In particular, we'll illustrate the method to approximate $\sqrt{m}$.

\subsection{Newton's Method}

A classic numeric problem is: Given a function
$f:\set{R}\mapsto\set{R}$, compute roots (zeros) $x$ such that $f(x)=0$.
In Western society, \Newton's name is attached to one method for
computing roots.
But the roots of the method can be traced to ancient mathematicians of \Mesopotomia,
the region of modern day Iraq and Iran.
Clay tablets from Mesopotamia, dated between 1800 B.C.~to~1600 B.C.,
showing   how to approximate \SQRootTwo and perform other arithmetic operations
\hlink{have been found}{http://www.math.ubc.ca/people/faculty/cass/Euclid/ybc/ybc.html}.
Algorithms that compute numbers like $\sqrt{2}$ or $\pi$ do not
terminate:  There is no finite positional notation for these numbers.

Newton is credited with generalizing the nascent idea to a broad class
of functions. An example function is $f(m)=x^{2}-m$, which has root $x=\sqrt{m}$.
The idea is: Let $x$ be a zero of function $f$.
Use \TaylorsTheorem.  Solve for $x$, and discard the second order
error.
\begin{align*}
  0 & = f(x) \approx f(x_{k}) + f'(x_{k})(x-x_{k}) + \frac{f''(\xi)}{2}(x-x_{k})^{2}\\
  x & = x_{k} - \frac{f(x_{k})}{f'(x_{k})} -\frac{f''(\xi)}{2f'(x_{k})}(x-x_{k})^{2}\\
  x & \approx x_{k} - \frac{f(x_{k})}{f'(x_{k})}
\end{align*}
For the function $f(x)=x^{2}-2$ the approximation reduces to 
\[  x \approx x_{k} - \frac{x_{k}^{2}-2}{2x_{k}}= \frac{x_{k}+2/x_{k}}{2}\]
The approximation is the average of $x_{k}>0$ and $2/x_{k}$.
If $x_{k}$ is too small: $x_{k}^{2}-2<0$, then 
$2/x_{k}$ is too large:
\[\bcond{(\bcond{(0 < 2 - x_{k}^{2})}{(0 < 2/x_{k}^{2} - 1)})}{((2/x_{k})^{2}-2>0)}\]
The average of $x_{k}$ and $2/x_{k}$ is a better approximation.

\begin{itemize}
\item Start with $h_{0}=1$ as an initial approximation to $\SQRootTwo$
\item Clearly $1$ is too small, as the \Babylonians could easily measure
\item But, if $1=\SQRootTwo$, then $1\cdot1=\SQRootTwo\cdot\SQRootTwo=2$
  and so $2/1=\SQRootTwo$
\item As it is, $2/1$ is too large
\item The average of the under estimate $1$ and the over estimate $2/1$
  provides a better approximation to $\SQRootTwo$, call this
  \[h_{1}=\frac{1}{2}\left(h_{0}+\frac{2}{h_{0}}\right)=\frac{1}{2}\left(1+\frac{2}{1}\right) = \frac{3}{2}\] 
\item But $h_{1}=3/2$ is too large, as the \Babylonians could measure
\item But, if $3/2$ were the exact square root, then $2/(3/2)=4/3$ would
  equal $\SQRootTwo$
\item As it is, $4/3\approx1.333\cdots$ is too small
\item The average of the over estimate $3/2$ and the under estimate $4/3$
  will provide a better approximation
  \[h_{2}=\frac{1}{2}\left(h_{1}+\frac{2}{h_{1}}\right)=\frac{1}{2}\left(\frac{3}{2}+\frac{4}{3}\right) = \frac{17}{12}\approx1.41166\cdots\]
\item But $h_{2}=17/12$ is too small
\item The \Babylonians carried out this iteration more times
  computing the $\SQRootTwo$ accurately to at least 9 decimal places
\item That is, they next computed the average of $h_{2}=17/12$ and $2/h_{2}=24/17$
  \begin{align*}
    h_{3} & =\frac{1}{2}\left(h_{2}+\frac{2}{h_{2}}\right)\\
    & =\frac{1}{2}\left(\frac{17}{12}+\frac{24}{17}\right)\\
    & =\frac{1}{2}\left(\frac{289 +288}{17\times12}\right) \\
    & \approx1.41421568628\cdots
  \end{align*}
\end{itemize}
To generalize the $\sqrt{2}$ method, 
pretend you want to compute $\sqrt{m}$. 
This is equivalent to computing a solution to the equation
\[x^{2}-m=0\]
Consider the recurrence
equation
\begin{equation}
  \label{newton:eq:1}
  x_{k} = \frac{x_{k-1} + m/x_{k-1}}{2},\quad k\geq1
\end{equation}
Given an initial value $x_{0}$, equation~\ref{newton:eq:1} can be used to
generate a sequence
\[\sequence{x_{0},\,x_{1},\,x_{2},\ldots}.\]
If you pretend that $x_{k}$ converges to $x$ as $k$ goes to infinity,
that is,
\[\lim_{k\rightarrow\infty}x_{k}=x\]
Then $x$ satisfies the equations
\begin{align*}
  x  & = \frac{x + m/x}{2} \\
  2x & = x + m/x \\
  x  & = m/x \\
  x^{2}  & = m \\
  x  & = \sqrt{m}
\end{align*}
%\todo[noline]{What happens when $m$ is negative?}

\subsection{Implementing Newton's method in Haskell}
\index{Listings!Haskell!Newton's Square Root Recurrence}

A first step in implementing Newton's method for computing $\sqrt{m}$ is to
define the function 
that maps $m$ and $x_{k-1}$ to the next value $x_{k}$.

\nwfilename{Newton.nw}\nwbegincode{1}\sublabel{NWNew9-New6-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWNew9-New6-1}}}\moddef{Newton~{\nwtagstyle{}\subpageref{NWNew9-New6-1}}}\endmoddef
module Newton ( newton ) where
\LA{}Apply the recursion~{\nwtagstyle{}\subpageref{NWNew9-AppJ-1}}\RA{}
\LA{}Repeatedly apply $f$ to $x$~{\nwtagstyle{}\subpageref{NWNew9-RepR-1}}\RA{}
\LA{}Test if successive values meet a relative tolerance~{\nwtagstyle{}\subpageref{NWNew9-Tesp-1}}\RA{}
\LA{}Newton's square root~{\nwtagstyle{}\subpageref{NWNew9-NewK-1}}\RA{}
\nwnotused{Newton}\nwendcode{}\nwbegindocs{2}\nwdocspar
The recursion is 
\[x_{k+1} = x_{k} - \frac{f(x_{k})}{f'(x_{k})}\quad\text{where $f(x)=x^{2}-2$}\]
\nwenddocs{}\nwbegincode{3}\sublabel{NWNew9-AppJ-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWNew9-AppJ-1}}}\moddef{Apply the recursion~{\nwtagstyle{}\subpageref{NWNew9-AppJ-1}}}\endmoddef
next :: Double -> Double -> Double
next m 0 = error "Division by zero"
next m x = (x + m/x)/2
\nwused{\\{NWNew9-New6-1}}\nwendcode{}\nwbegindocs{4}\nwdocspar

We want to repeatedly apply \lstinline{next} to some initial value and
generate a list of \lstinline{Double}s.  Let's define
\lstinline{repeatedly} to be a function that applies a function
\lstinline{f :: Double -> Double} to itself repeatedly.  An initial
value (seed) \lstinline{a} for {\tt{}f} starts the iteration, generating
an infinite list.

\index{Listings!Haskell!Apply a Function Repeatedly}
\nwenddocs{}\nwbegincode{5}\sublabel{NWNew9-RepR-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWNew9-RepR-1}}}\moddef{Repeatedly apply $f$ to $x$~{\nwtagstyle{}\subpageref{NWNew9-RepR-1}}}\endmoddef
repeatedly :: (Double -> Double) -> Double -> [Double] 
repeatedly f a = a : repeatedly f (f a)
\nwused{\\{NWNew9-New6-1}}\nwendcode{}\nwbegindocs{6}\nwdocspar

Although \lstinline{repeatedly} does not terminate,
it can be terminated once a computed value is \emph{close enough}.
A common way to do this is to define a tolerance %\tau$,
usually the \machineepsilon and declare that the 
last computed approximation is \emph{good enough} once it 
and the previous approximation are within the tolerance.

The \emph{absolute difference} $\abs{x_{k}-x_{k-1}}$ between successive
iterates is a measure of closeness.
One way to terminate Newton's iteration is to stop when the absolute
difference is within the tolerance.
\[\abs{x_{k}-x_{k-1}}\leq\tau\]
Computer arithmetic on floating point numbers is not exact.  The
absolute difference can be small because the numbers $x_{k}$ and
$x_{k-1}$ themselves are small. The absolute difference may never be
small because the numbers themselves are large.


Instead of computing until the difference of successive approximations
approaches $0$, it is often better to compute until the ratio of
successive approximations approach $1$.  This measure of closeness is the \emph{relative
  difference} $\abs{x_{k-1}/x_{k}-1}$.
Some define the relative error as
\[re_{k} = \frac{x_{k} - x}{x}\]
where $x$ root being sought.  But this requires knowledge of $x$ to compute.

Newton's method terminates when 
\[\abs{\frac{x_{k-1}}{x_{k}}-1}\leq\tau\]
The {\tt{}relative} function maps a tolerance $\tau$ and a sequence to
the first value in the sequence where the relative error is within
tolerance.
%\begin{mylisting}{Convergence of Relative Error}{rel-error}
\index{Listings!Haskell!Relative Error Convergence Test}
\nwenddocs{}\nwbegincode{7}\sublabel{NWNew9-Tesp-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWNew9-Tesp-1}}}\moddef{Test if successive values meet a relative tolerance~{\nwtagstyle{}\subpageref{NWNew9-Tesp-1}}}\endmoddef
relative :: Double -> [Double] -> Double
relative tau (a:b:rest) 
   | abs (a/b-1) <= tau = b
   | otherwise          = relative tau (b:rest)
\nwused{\\{NWNew9-New6-1}}\nwendcode{}\nwbegindocs{8}\nwdocspar
%\end{mylisting}
\end{fullwidth}

Now we can express Newton's method to compute the square root of
$m$ to within a relative error tolerance $\tau$ starting with an initial guess
$x_{0}$ as the function {\tt{}newton}.
\todo[noline]{Show how the code can be changed to compute other
  values. Test this on other functions, for instance $f(x)=x^{3}-2$ to
  compute $\sqrt[3]{2}$.}
%\begin{mylisting}{Newton's Square Root Method}{newton-sqroot}
\index{Listings!Haskell!Newton's Square Root Method}
\nwenddocs{}\nwbegincode{9}\sublabel{NWNew9-NewK-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWNew9-NewK-1}}}\moddef{Newton's square root~{\nwtagstyle{}\subpageref{NWNew9-NewK-1}}}\endmoddef
newton x0 tau m = relative tau (repeatedly (next m) x0)
\nwused{\\{NWNew9-New6-1}}\nwendcode{}

\nwixlogsorted{c}{{Apply the recursion}{NWNew9-AppJ-1}{\nwixu{NWNew9-New6-1}\nwixd{NWNew9-AppJ-1}}}%
\nwixlogsorted{c}{{Newton}{NWNew9-New6-1}{\nwixd{NWNew9-New6-1}}}%
\nwixlogsorted{c}{{Newton's square root}{NWNew9-NewK-1}{\nwixu{NWNew9-New6-1}\nwixd{NWNew9-NewK-1}}}%
\nwixlogsorted{c}{{Repeatedly apply $f$ to $x$}{NWNew9-RepR-1}{\nwixu{NWNew9-New6-1}\nwixd{NWNew9-RepR-1}}}%
\nwixlogsorted{c}{{Test if successive values meet a relative tolerance}{NWNew9-Tesp-1}{\nwixu{NWNew9-New6-1}\nwixd{NWNew9-Tesp-1}}}%
\nwbegindocs{10} 
%\end{mylisting}

\subsection{Convergence of Newton's Method}

\begin{fullwidth}

The number of times the \lstinline{(next m)} function is evaluated 
measures the time complexity of the \lstinline{mysqrt} algorithm.
It is not obvious what this number is.  
What can be shown is that Newton's method converges quadratically,
under certain assumptions that are often \true.

There is a more elegant, more general proof of the
quadratic convergence of Newton's method.  But this is beyond the
scope of these notes.



Consider the function $f(x)=x^{2}-m$.  
Using \TaylorsTheorem, you can derive the equation
\[x^{2}-m = (x_{k-1}^{2}-m) + 2x_{k-1}(x-x_{k-1})+(x-x_{k-1})^{2}\]
Pretend that $x=\sqrt{m}$ so that both sides of the above equation are
zero. Divide by $2x_{k-1}$ to get
\[0  = (x_{k-1}^{2}-m)/2x_{k-1} + (x-x_{k-1})+(x-x_{k-1})^{2}/2x_{k-1}\]
Notice that
\[(x_{k-1}^{2}-m)/2x_{k-1} - x_{k-1} = -(x_{k-1}^{2}+m)/2x_{k-1} = -x_{k}\]
Therefore,
\begin{align*}
  x_{k}-x & = (x_{k-1}-x)^{2}/2x_{k-1}\quad\text{or}\\
  e_{k} & = \frac{e_{k-1}^{2}}{2x_{k-1}} 
\end{align*}
That is, the absolute error $x_{k}-x$ at step $k$ is proportional to the square
of the error at step $k-1$.  When the error is less than $1$, the
number of correct digits doubles with each iteration.

In the general case, assume that function $f$ has a continuous second
derivative.  Assume $x$ is a root of $f$, that is $f(x)=0$.
By Taylor's theorem 
  \begin{align*}
  0 & = f(x) \approx f(x_{k}) + f'(x_{k})(x-x_{k}) + \frac{f''(\xi)}{2}(x-x_{k})^{2}\\
  x & = x_{k} - \frac{f(x_{k})}{f'(x_{k})} -\frac{f''(\xi)}{2f'(x_{k})}(x-x_{k})^{2}\\
  x & \approx x_{k} - \frac{f(x_{k})}{f'(x_{k})}
  \end{align*}
  
\end{fullwidth}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: t
%%% End: 
\nwenddocs{}
