\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{The bisection method}\\
  \Large{MATH2070: Numerical Methods in Scientific Computing I}
\end{center}
Location: http://people.sc.fsu.edu/$\sim$jburkardt/classes/math2070\_2019/nonlinear\_bisection/nonlinear\_bisection.pdf
\vskip 0.1in
\hrule
\vskip 0.1in
\begin{center}
  \includegraphics[scale=0.18]{apple_split.png}
  \vskip 0.1in
  {\it{Why does a function cross the axis? To get to the other sign!}}
\end{center}

Mathematics provides methods for finding solutions to some algebraic
equations like $x^2-2x-15=0$.  However, for most equations with any kind of
complication, there is no way to come to an exact solution for $x$.
The bisection method provides a computational path to solving a nonlinear
equation.
\vskip 0.1in
\begin{tcolorbox}[colback=red!5,title=The bisection method]
{\it{Given a nonlinear equation:
\begin{packed_itemize}
  \item{rewrite it as $f(x)=0$}
  \item{find points $a$ and $b$ at which $f(x)$ takes opposite signs;}
  \item{shrink this interval until you are close enough to a solution.}
\end{packed_itemize}
}}
\end{tcolorbox}

\section{Example: Kepler's equation}

Kepler's equation comes from an astronomical problem.  It 
relates an important quantity, $E$, the eccentric anomaly of an orbit, to
two easily measured items, the mean anomaly $M$ and the orbital 
eccentricity $e$:
\begin{displaymath}
  M = E - e \sin(E)
\end{displaymath}
\begin{center}
  \includegraphics[scale=0.40]{kepler_equation.png}
  \vskip 0.1in
  \it{$E$ is the ``height'' difference between an ideal and elliptic orbit. }
\end{center}
Although we may know specific values of $M$ and $e$, it is not possible 
to rewrite this equation to determine the exact value of $E$ except for
a few special cases.  Suppose that we have determined that $M = 5$ and $e = 2$.  To determine the value
of $E$, we need to solve the equation:
\begin{displaymath}
  5 = E - 2 \sin(E)
\end{displaymath}
We rewrite this equation in functional form:
\begin{displaymath}
  f(E) = 5 - E + 2 \sin(E) = 0
\end{displaymath}
Assuming that $E^*$ is an exact solution of this equation, we have three 
common measures of error:
\begin{packed_itemize}
  \item{$|E^*-E|$ is the (absolute) error;}
  \item{$\frac{|E^*-E|}{|E^*|}$ is the relative error;}
  \item{$|f(E)|$ is the (absolute) residual error;}
\end{packed_itemize}
\vskip 0.1in
To begin our investigation, let's write a MATLAB function to evaluate the residual:
\begin{lstlisting}[caption=kepler.m]
function value = kepler ( x )
  value = 5.0 - x + 2.0 * sin ( x );
  return
end
\end{lstlisting}
We use plotting to search for a region in which the function
crosses the axis:
\begin{center}
  \includegraphics[scale=0.40]{kepler.png}
  \vskip 0.1in
\end{center}
By evaluation, we see that $\text{kepler}(0)=5$ and $\text{kepler}(10)=-6.0880$.  
As $E$ moves from 0 to 10, $f(E)$ moves from positive to negative. 
Assuming $f(E)$ is continuous, there must be a solution $E^*$ 
somewhere in $[0,10]$.  Looking only at the change in sign, our best guess might 
be that the solution is near $E=5$.
\vskip 0.1in
In fact, $\text{kepler}(5) = -1.9178$,  which suggests we can find the root 
in $[0,5]$.  And if we keep looking at the midpoint of our current interval,
we can repeatedly cut our uncertainty in half, and get as close to the root
as we want (up to the limits of the machine accuracy).

\section{Pseudocode \#1: Bisection}

Assume we know $xn$ and $xp$ so that $f(xn)<0$ and
$f(xp)>0$.  Our bisection algorithm might look like:
\begin{lstlisting}[caption = bisection1\_pseudocode.txt]
Input: function f, negative and positive arguments x- and x+

Loop a few times

  Set x to the average of x- and x+
  
  If f(x) is negative, replace x- by x
  else replace x+ by x

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

\section{MATLAB implementation \#1:}

Our simple bisection algorithm needs to stop at some point, so
let's just let it run 10 steps:
\begin{lstlisting}[caption=bisection1.m]
function [ xn, xp ] = bisection1 ( f, xn, xp )

  for it = 1 : 10

    x = ( xn + xp ) / 2.0;

    if ( f ( x ) < 0.0 )
      xn = x;
    else
      xp = x;
    end
    
  end

  return
end
\end{lstlisting}

\section{Pseudocode \#2: Bisection}

In class, we talked about sensible ways to tell an iteration to
continue, or to stop with success or failure.  This means we need
to supply an iteration limit, and tolerances on $x$ and $f(x)$.
Our second pseudocode might be:
\begin{lstlisting}[caption = bisection2\_pseudocode.txt]

Input: function f, negative and positive arguments x- and x+, xtol, ftol, itmax

Begin loop

  Set x to the average of x- and x+
  
  If f(x) is negative, replace x- by x
  else replace x+ by x
  
  new = | x+ - x- |

  if new is less than xtol and |f| is less than ftol, success
  else if it > itmax, failure
  else repeat

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

\section{MATLAB implementation \#2:}

Our revised code might be:
\begin{lstlisting}[caption=bisection2.m]
function [ xn, xp, it ] = bisection2 ( f, xn, xp, xtol, ftol, itmax )

  it = 0;

  while ( true )

    it = it + 1;

    x = ( xn + xp ) / 2.0;

    if ( f ( x ) < 0.0 )
      xn = x;
    else
      xp = x;
    end

    new = abs ( xn - xp )

    if ( ...
      ( new <= xtol ) && ...
      ( abs ( f ( x ) ) <= ftol ) ...
    )
      return
    end

    if ( itmax <= it )
      return
    end

  end

  return
end
\end{lstlisting}

\section{Bisection for the Kepler equation}

Test bisection code \#2 on our Kepler equation:

\begin{lstlisting}[caption =  kepler\_bisection.m]
xn = 10.0;
xp = 0.0;
xtol = 0.000001;
ftol = 0.000001;
itmax = 50;

[ xn, xp, it ] = bisection2 ( @(x)kepler(x), xn, xp, xtol, ftol, itmax );
\end{lstlisting}

Print statements in the script (not shown) report the following results:

\begin{lstlisting}[caption = kepler\_bisection.txt]
  After 24 iterations:
  F(3.79084) = -1.52126e-06 (negative)
  F(3.79083) = 2.43377e-08 (positive)

  it <= itmax?  true
  |xn-xp| <= xtol?  true
  min(|f(xn)|,|f(xp)|) <= ftol?  true
\end{lstlisting}

To apply the bisection code to another problem, write a 
function file to evaluate $f(x)$, and modify the test script by
specifying appropriate points $xn$ and $xp$.  You might also want 
to adjust the tolerances or the iteration limit.

\section{Exercise: A trigonometric function}

Consider the equation$\cos(x) = x$
and suppose we want to find a solution $x$.  Use the bisection method
to approximate such a solution.  To do this:
\begin{packed_enumerate}
  \item{create a function file {\it{trig.m}} that evaluates the function.}
  \item{create a script file like {\it{trig\_bisection.m}} that calls
  {\tt{bisection2()}} to find a zero.}
  \item{run your script;}
  \item{report the number of iterations required;}
  \item{report as your solution $x$ the average of the two endpoints;}
  \item{report the value of $f(x)$;}
\end{packed_enumerate}

\section{MATLAB implementation \#3: ALPHA, the update ratio:}

In class, the quantity $\alpha$ or ``alpha'' was defined as the
ratio $\alpha = \frac{\text{this x update}}{\text{previous x update}}$.
For the bisection method, we can think of $\alpha$ as the ratio of
the new interval to the old one, so that for bisection, $\alpha$ is 
always $\frac{1}{2}$.  Nonetheless, here is how we could go through
the motions of computing $\alpha$:
\begin{lstlisting}[caption=bisection3.m]
function [ xn, xp, it ] = bisection3 ( f, xn, xp, xtol, ftol, itmax )

  it = 0;

  while ( true )

    old = abs ( xp - xn );
    
    it = it + 1;

    x = ( xn + xp ) / 2.0;
    
    if ( f ( x ) < 0.0 )
      xn = x;
    else
      xp = x;
    end

    new = abs ( xp - xn );
    alpha = new / old;
    
    if ( ...
      ( abs ( xn - xp ) <= xtol ) && ...
      ( abs ( f ( x ) ) <= ftol ) ...
    )
      return
    end

    if ( itmax <= it )
      return
    end

  end

  return
end
\end{lstlisting}
\vskip 0.1in

\section{Fixed point iteration}

In class, we saw a general method to solve nonlinear
equations, called {\it{fixed point iteration}}.  To find $x^*$ that satisfies
$f(x)=0$, we repeatedly update $x \leftarrow g(x)$, 
where $g(x)$ is chosen so that $x^*=g(x^*)$.
\vskip 0.1in
We found a root of $f(x) = 3 * x - e^x = 0$ by using the fixed point 
function $g(x) = e^x/3$.  Pseudocode for this computation would be:
\begin{lstlisting}[caption = fixed\_point\_pseudocode.txt]
itmax <-- ?
xtol <-- ?
ftol <-- ?

it <-- 0
xold <-- 0
x <-- initial value   (1.0 for this case)
old <-- 0
new <-- 0

Loop 

  it <-- it + 1
  xold <-- x
  x <-- g(x)          (x = exp(x)/3; for this case)
  old <-- new
  new <-- abs ( x - xold )

  if old is not 0, alpha <-- new / old, print alpha
  if new <= xtol and |f(x)| <= ftol, break from loop with success
  if it > itmax, break from loop with failure

End loop
  
Print x, f(x)

function value = f ( x )
  value <-- ?;  ( value = 3 * x - exp ( x ); for this case )
  return
end
\end{lstlisting}

\section{Assignment \#3: Implement and test a fixed point iteration}

Use the  pseudocode above as a guide, write a script called
{\it{hw3.m}} which uses fixed point iteration to find a root of
the function $f(x)=x^2-5$.  Use the fixed point function 
$g(x)=1+x-x^2/5$.
\vskip 0.1in
{\bf{Turn in:}} your file {\it{hw3.m}} by Friday, September 13.

\end{document}
