\documentclass{article}

\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{graphics}
\usepackage{listings}
\usepackage{tcolorbox}
\usepackage{verbatim}

\DeclareMathOperator*{\argmin}{arg\,min}

\marginparwidth 0pt
\oddsidemargin  0pt
\evensidemargin  0pt
\marginparsep 0pt
\topmargin   0pt
\textwidth   6.5in
\textheight  8.5in
\setlength{\parindent}{0em}

\newenvironment{packed_itemize}{
\begin{itemize}
  \setlength{\itemsep}{0pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{itemize}}

\newenvironment{packed_enumerate}{
\begin{enumerate}
  \setlength{\itemsep}{1pt}
  \setlength{\parskip}{0pt}
  \setlength{\parsep}{0pt}
}{\end{enumerate}}

\definecolor{codeback}{RGB}{235,255,255}

\lstset{ %
language=Matlab,                % choose the language of the code
basicstyle=\footnotesize,       % the size of the fonts that are used for the code
numbers=left,                   % where to put the line-numbers
numberstyle=\footnotesize,      % the size of the fonts that are used for the line-numbers
stepnumber=1,                   % the step between two line-numbers. If it is 1 each line will be numbered
numbersep=10pt,                  % how far the line-numbers are from the code
backgroundcolor=\color{codeback},  % choose the background color. You must add \usepackage{color}
showspaces=false,               % show spaces adding particular underscores
showstringspaces=false,         % underline spaces within strings
showtabs=false,                 % show tabs within strings adding particular underscores
frame=single,           % adds a frame around the code
tabsize=2,          % sets default tabsize to 2 spaces
captionpos=b,           % sets the caption-position to bottom
breaklines=true,        % sets automatic line breaking
breakatwhitespace=false,    % sets if automatic breaks should only happen at whitespace
escapeinside={(*@}{@*)}          % if you want to add a comment within your code
}

\begin{document}

\begin{center}
  \Huge{Nonlinear equations:}\\
  \Huge{Newton's method}\\
  \Large{MATH2070: Numerical Methods in Scientific Computing I}
\end{center}
Location: http://people.sc.fsu.edu/$\sim$jburkardt/classes/math2070\_2019/nonlinear\_newton/nonlinear\_newton.pdf
\vskip 0.1in
\hrule
\vskip 0.1in
\begin{center}
  \includegraphics[scale=0.30]{ski_slope.png}
  \vskip 0.1in
  {\it{Even in zero visibility, a skier knows which way is downhill!}}
\end{center}

\begin{tcolorbox}[colback=red!5,title=Rootfinding (Newton version)]
{\it{Given $f(x)$, with derivative $f'(x)$, and
a starting point $x$, estimate a root $x^*$
for which $f(x^*)=0$?}}
\end{tcolorbox}

\section{Newton's linear model for f(x) }

Suppose that we are looking for a root $x^*$ of $f(x)=0$, and that 
$f(x)$ is at least twice continuously differentiable.  Then, at a
point $x$ sufficiently close to $x^*$, we can write:
\begin{displaymath}
f(x*) = f(x) + f'(x) * ( x^* - x ) + O(x^*-x)^2
\end{displaymath}
Noting that $f(x^*)=0$, and dropping the $O(x^*-x)^2$ terms, we arrive
at Newton's estimate for the location of the root:
\begin{displaymath}
x* \approx x - \frac{f(x)} {f'(x)}
\end{displaymath}
\vskip 0.1in
Thus, to estimate the root of a function, we try to find a starting
point that we hope is already close, and we must be prepared to
evaluate the derivative of the function as part of the iteration.
\vskip 0.1in
We hope that, once the iteration is close to the root, the norm 
of the function will decrease on every step, and the size of the
step $|x^i-x^{i-1}|$ will decrease quadratically, as it becomes
similar to the size of the actual (but unknowable) error $|x^*-x^{i-1}|$.

\section{Newton's method pseudocode}

The pseudocode for Newton's method is very similar to that for
previous algorithms.  However, we only pass in a single starting
value {\tt{x}}, and along with the function {\tt{f(x)}} we now
must include the derivative function {\tt{fp(x)}}.

\begin{lstlisting}[caption = Pseudocode for Newton's method]
## newton_pseudocode.txt

function newton ( )

Input: f(), (*@ \textcolor{red}{fp(),} @*) x, xtol, ftol, itmax

it = 0
dx = 0
alpha = 0

Begin loop

  it = it + 1

  Set xold to x, and set (*@ \textcolor{red}{x to x - f(x) / fp(x)} @*)

  Set dx_old to dx, and set dx = | x - xold |

  Set alpha_old to alpha
  if dx_old is not 0, set alpha to dx / dx_old 
    if alpha_old is not 0, set r = log ( alpha ) / log ( alpha_old )

  if x - xold is less than xtol and |f(x)| is less than ftol, success
  if it > itmax, failure

End loop
  
Output: values of updated x and it
\end{lstlisting}

\section{Example: Kepler's Equation}

Our example of a Kepler function has the form:
\begin{displaymath}
  f(x) = 5 - x + 2 * \sin ( x )
\end{displaymath}
so the corresponding derivative function would be:
\begin{displaymath}
  f'(x) = -1 + 2 * \cos ( x )
\end{displaymath}
To use Newton's method, in addition to our file {\it{kepler.m}}, 
we would need to supply an additional file {\it{kepler\_fp.m}} of
the form:
\begin{lstlisting}[caption = kepler\_fp.m]
function value = kepler_fp ( x )

  value = - 1.0 + 2.0 * cos ( x );

  return
end
\end{lstlisting}
\vskip 0.1in
Let us apply our Newton code to this problem.  We can call our
MATLAB Newton solver with this script:
\begin{lstlisting}[caption = Newton's method applied to the Kepler function]
x = 3.0;
xtol = 0.00000001;
ftol = 0.00000001;
itmax = 50;

[ x, it ] = newton ( @(x)kepler(x), @(x)kepler_fp(x), x, xtol, ftol, itmax );
\end{lstlisting}
The algorithm returns after just four steps.  Here are the intermediate
results:
\begin{center}
\begin{tabular}{rll}
it&      x&                 f(x)\\
\hline
0&     3.000000000000000&   2.28224\\
1&     3.765856211130886&   0.065144\\
2&     3.790693912837725&   0.00036471\\
3&     3.790834550863825&   1.19566e-08\\
4&     3.790834555474780&  -2.22045e-16\\
\end{tabular}
\end{center}
While the values $f(x)$ seem to be showing quadratic convergence,
the true quadratic behavior is in the successive values of $x$,
in which the number of correct digits seems to roughly double
on each step.
  
\section{Exercise: The jump function}

As an in-class exercise, let's see what we need to do to apply
Newton's method to a new example, the {\it{jump}} function,
whose plot looks a little like a ski jump.  
The function is defined by:
\begin{displaymath}
  f(x)=e^{\cos(x)}-x/2+2
\end{displaymath}
\vskip 0.1in
Do the following:
\begin{packed_enumerate}
  \item{Copy the file {\it{newton.m}} from the website;}
  \item{Create a function file {\it{jump.m}};}
  \item{Create a derivative file {\it{jump\_fp.m}};}
  \item{Plot $f(x)$ between 0 and 10.  Use this plot to pick
  a good starting point $x$.  Choose $x$ to be an integer;}
  \item{Write a script, or interactively issue the commands,
  which specify the input to the {\it{newton()}} function;}
  \item{Report your starting point, your solution, and the number of iterations required.}
\end{packed_enumerate}
\vskip 0.1in
As with the secant method, our best results occur if the starting
point is in the region near the root where the function seems to
be behaving linearly.  If the starting point is further away, then
Newton's method can be very slow, or wander aimless, or even diverge
to infinity.

\section{Example: Convergence rate for the Lambert function}

Our Lambert function has the form
\begin{displaymath}
f(x) = x * e^x - 742.0657955128830
\end{displaymath}
For Newton's method, we also need the derivative:
\begin{displaymath}
f'(x) = ( x + 1 ) * e^x
\end{displaymath}
We can use a starting value of $x=3$ and call the {\tt{newton()}} 
function to seek a root.  But now, we are interesting in observing
the convergence rate of the algorithm.
\vskip 0.1in
As we have done for previous algorithms, we can estimate the
order of convergence $r$ as we come near the root  by comparing
the differences between successive iterates.  Newton's method should
in general exhibit quadratic convergence, so the value of $r$
should tend to 2.  
\vskip 0.1in
If we turn on the print statements inside {\tt{newton()}} which
monitor the convergene rate estimate, here is what we see for
a calculation with the Lambert function:
\begin{center}
\begin{tabular}{rlll}
  it & $\alpha$ & $\log(\alpha)$ & r \\
  \hline
  2 &  0.108328 &    -2.2226 & \\
  3&  0.991943&    -0.00808959& 0.00363971\\
  4&  0.988803&    -0.0112604& 1.39197\\
  5&  0.982003&    -0.0181604& 1.61276\\
  6&  0.965831&    -0.0347668& 1.91442\\
  7&  0.925889&    -0.0770008& 2.21478\\
  8&  0.83009&     -0.186221& 2.41843\\
  9&  0.630036&    -0.461978& 2.48081\\
 10&  0.331472&    -1.10421& 2.39018\\
 11&  0.0889238&   -2.41998& 2.19159\\
 12&  0.00718688&  -4.9355& 2.03948\\
 13&  5.11763e-05& -9.88023& 2.00187\\
\end{tabular}
\end{center}
Looking at the values of $\alpha$, it's not so easy to see quadratic
convergence, but the $\log(\alpha)$ column makes it easier to see.

\section{Exercise: The multi function}

Consider the problem of finding a root of another new example,
called the ``multi'' function.  This function has the form
\begin{displaymath}
  f(x) = x^4 - 6x^2 + 8x -3
\end{displaymath}
with derivative
\begin{displaymath}
  f'(x) = 4x^3 - 12x + 8
\end{displaymath}
\begin{packed_itemize}
  \item{Plot the function over [-3.5,+3.0], and notice there seem to be
  two places where $f(x)=0$;}
  \item{Using the starting point $x=3$, call the {\tt{newton()}} function
  to get a root.  Note the number of iterations seems high;}
  \item{Turn on the print statement in {\tt{newton()}} which reports
  the convergence rate $r$.  What seems to be different for this problem?}
  \item{The actual value of the root is an integer, so if {\tt{newton()}}
  didn't quite get there, round the result to an integer.  Verify that
  $f(x)$ is exactly 0 at this value of $x$.}
  \item{Now verify that it is also
  true that $f'(x)=0$.  Look at the plot again to see that this
  is so.  When a function and its derivative are both zero at a value
  $x$, what does this say about $x$?}
  \item{What was said in class about the behavior of Newton's method
  in this case?}
\end{packed_itemize}

\section{No assignment for this lab!}

\end{document}
