\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{Differentiation:}\\
  \Huge{The Finite Difference Method}\\
  \Large{MATH2070: Numerical Methods in Scientific Computing I}
\end{center}
Location: http://people.sc.fsu.edu/$\sim$jburkardt/classes/math2070\_2019/differentiation\_difference/differentiation\_difference.pdf
\vskip 0.1in
\hrule
\vskip 0.1in

\begin{center}
  \includegraphics[scale=0.50]{differences}
  \vskip 0.1in
  {\it{Difference estimates for the derivative}}
\end{center}

\vskip 0.1in
\begin{tcolorbox}[colback=red!5,title=Estimating Derivatives]
{\it{
  Given the formula for a function $f(x)$:
  \begin{packed_itemize}
    \item{how can we estimate derivatives, in particular, $f'(x)$?}
    \item{how does the finite difference stepsize affect the error?}
    \item{how does roundoff error limit the achievable accuracy?}
  \end{packed_itemize}
}}
\end{tcolorbox}
\vskip 0.1in

\section{The forward difference estimate to $f'(x)$:}

Suppose we have a formula, or a computational procedure, that provides us with
the value $f(x)$ of a function, but now we need to estimate the derivative of the
function, and calculus is not an option.  The very definition of the derivative
says it is the limit of the ratio of differences in function values to differences
in argument values.  Therefore, if $h$ is ``small enough'', we can expect
a good estimate for $f'(x)$ using the forward difference, symbolized by 
{\bf{fd(f,x,h)}}:
\begin{displaymath}
  f'(x) \approx \text{fd}(f,x,h) \equiv \frac{f(x+h)-f(x)}{h}
\end{displaymath}
It's easy to create a corresponding MATLAB function that will produce the 
forward difference estimate:
\begin{lstlisting}
function value = fd ( f, x, h )

  value = ( f ( x + h ) - f ( x ) ) / h;

  return
end
\end{lstlisting}
Notice that the first argument is not a number, but the function.  When using
{\tt{fd()}}, then, we can't simply pass a formula or function name; we have
to identify this argument as a function, and specify its argument.  Here
are a few examples:
\begin{lstlisting}
  value = fd ( @(x)tan(x), x, h );
  value = fd ( @(x)log(x+1), x, h );
  value = fd ( @(x)x^2+3*x+4, x, h );   % A formula
\end{lstlisting}
Practice this command by setting {\tt{x=1}} and {\tt{h=0.1}} and then trying
each of the three examples above.  You should see answers of 4.0735, 0.4879, and
5.1.
\vskip 0.1in
Notice that {\tt{fd()}} can also be called with a vector {\tt{x}}.  That means
we can make plots comparing the exact and estimated derivatives.  Try the
following commands:
\begin{lstlisting}
  x = linspace ( 0.0, 2.0, 101 )   % Remember this command?
  h = 0.1;
  y = fd (@(x)log(x+1), x, h );
  y2 = 1 ./ ( x + 1 );
  plot ( x, y, 'r-', x, y2, 'b-' );  % Draw red and blue graphs together.
\end{lstlisting}

What ``breaks'' if you repeat this computation with the command:
\begin{lstlisting}
value = fd ( @(x)x^2+3*x+4, x, h );
\end{lstlisting}
Can we fix the problem?

\section{Discretization and roundoff errors}

In exact arithmetic, the difference between the true and estimated values 
is known as the {\bf{discretization error}}.  We can estimate the order of the 
absolute value of the discretization error:
\begin{displaymath}
  |f'(x) - \text{fd}(f,x,h)| = O(h)
\end{displaymath}
\vskip 0.1in
This formula seems to promise that we can drive the finite difference estimate
towards the true derivative value simply by relentless decreasing $h$.  But this
is only half of the story!
\vskip 0.1in
As explained in class, as $h$ is decreased, we began to make increasingly
dangerous arithmetic errors.  The quantities $f(x+h)$ and $f(x)$ become
very close, so the number of significant digits in their difference decreases, and
so we compute something slightly different from the exact arithmetic forward difference.
This is known as {\bf{roundoff error}}.
As $h$ gets smaller, the discretization error decreases, but the roundoff error
rises.  This means that at some point, we can expect our total error to begin
to increase, and after that, our computed estimates are worthless.

\section{Compare the forward difference estimate to the exact derivative}

In class, the example of estimating the derivative of the tangent function
at $x=\frac{\pi}{4}$ was used.  We know that $tan'(\frac{\pi}{4})=2$.  We 
imagine that decreasing $h$ gives us better and better estimates.  Let's see
if this is true.
\vskip 0.1in
For convenience, let's modify our forward difference function so that it can 
compute the estimates for a vector of $h$ values:
\begin{lstlisting}
function value = fd ( f, x, h )

  value = ( f ( x + h ) - f ( x ) ) ./ h;

  return
end
\end{lstlisting}
It's a tiny difference, but it will make things much easier.  Now our input
{\tt{h}} can be a list of steps of decreasing size, and we can compute all
the finite difference estimates at once.
\begin{lstlisting}
  h = 0.5;
  alpha = 0.5;
  n = 50
  p = 0 : n - 1;              % p = [ 0, 1, 2, ..., 49 ]
  h_list = h * alpha .^ p;    % h_list = [ h, h*alpha, h*alpha^2, ..., h*alpha^(n-1)
  fd_list = fd ( @(x)tan(x), pi/4, h_list )
  e_list = abs ( 2.0 - fd_list );
  plot ( h_list, e_list );
\end{lstlisting}
Looking at the plot, everthing seems great.
\vskip 0.1in
But if you actually look at the  {\tt{fd\_list}} values, you can see that
something odd is happening.  At first the values seem to converge towards 2,
but then they pull away again.  (You may want to try the {\tt{format long}} 
statement to see more digits in these values.)  Normally, a plot reveals
information, but this time, it seems to be hiding things from us.  What is
going wrong?
\vskip 0.1in
If you think about it, all the bad stuff is happening at extremely small
values of $h$, so in order to see the disaster, we would need a microscope 
to be able to examine the lower left corner of the plot.  Unfortunately,
the MATLAB plotting algorithm doesn't actually plot these values because
they are so close together!
\vskip 0.1in
When numbers vary greatly in scale, one way to try to fit them into a plot
is to work with the logarithm.
\begin{lstlisting}
  h = 0.25; 
  alpha = 0.5;
  n = 50;
  p = 0 : n - 1;
  h_list = h * alpha .^ p;
  fd_list = fd ( @(x)tan(x), pi/4, h_list );
  e_list = abs ( 2.0 - fd_list );
  %* {\color{red}plot ( log ( h\_list ), log ( abs ( e\_list ) ) ); } %
\end{lstlisting}
Now the information in our data can be displayed on the plot in a way that
tells a story!  We can use the {\tt{min()}} command to get the value (1.5e-08) 
and the index (26) of the minimum entry in {\tt{e\_list}}.  
\begin{lstlisting}
  [ e_min, i_min ] = min ( e_list );
  fprintf ( '  Minimum error is at index %d\n', i_min );
  fprintf ( '  Minimum error is %g\n', e_min );
  fprintf ( '  Best h is %g\n', h_list(i_min) );
  fprintf ( '  Best estimate is %20.16g\n', fd_list(i_min) );
\end{lstlisting}
This means that 
our best estimate for the derivative was {\tt{fd\_list(26)=2.000000014901161}}.  
It means that the minimum useful 
value of {\tt{h}} is {\tt{h\_list(26)}} or 7.45e-09, which is about half of 
the square root of the machine precision {\tt{sqrt(eps)=1.49e-08}}.

\section{Assignment \#2}

A better estimate for $f'(x)$ uses a centered difference.  We can symbolize this by 
{\bf{center(f,x,h)}}:
\begin{displaymath}
  f'(x) \approx \text{center}(f,x,h) \equiv \frac{f(x+h)-f(x-h)}{2h}
\end{displaymath}
We can estimate the order of the 
absolute value of the discretization error:
\begin{displaymath}
  |f'(x) - \text{center}(f,x,h)| = O(h^2)
\end{displaymath}
\vskip 0.1in
In class, it was suggested that this improved accuracy meant that the centered
difference error would continue to decrease down to a smaller $h$ than we saw
for the forward difference.
\vskip 0.1in
Write a corresponding MATLAB function {\tt{center.m}} to compute the centered
difference estimate of the first derivative.  Repeat the exercise for
estimating the first derivative of $\tan(x)$ at $x=\frac{pi}{4}$, using 50 values
of $h$, starting at $h=0.25$ and dividing by half each time.  
Determine the value of $h$ at which you made the smallest error, given that the correct derivative
is 2.
\vskip 0.1in
Your program, which might be called {\it{center\_test.m}}, should
\begin{packed_enumerate}
  \item{set up 50 decreasing values of $h$ in 
  {\tt{h\_list}};}
  \item{use your {\tt{center()}} function to compute 
  the centered difference estimates
  {\tt{cd\_list}} for all the values in {\tt{h\_list}};}
  \item{determine the value of $h$ in 
  {\tt{h\_list}} at which you made the 
  smallest error;}
  \item{plot {\tt{h\_list}} versus {\tt{cd\_list}};}
\end{packed_enumerate}
Send me your program via email.  I don't need a copy of your plot.
The assignment is due by Friday evening, September 6.

\end{document}
