\documentclass[12pt]{article}
\textwidth 16cm
\textheight 23cm
\topmargin 0cm
\evensidemargin 0cm
\oddsidemargin 0cm
\parindent=0pt
\parskip=2pt
\begin{document}


\begin{center}
\Large \bf
Computational Methods and Data Analysis 2003\\
\large 
\end{center}
\vspace{0.4cm}
\noindent 
{\bf Ex.~2: Numerical derivatives}\\

Make a new directory CMDA/ex2 and open it for reading.

The first and second derivative of a function can be approximated numerically
as finite differences.  From Taylor's expansions
one can derive the forward difference formula
\begin{equation}
  f'(x) = \frac{f(x + h) - f(x)}{h} + {\cal O}(h)
\end{equation}

as well as  the backwards difference formula

\begin{equation}
f'(x) \approx \frac{f(x) - f(x - h)}{h}
\end{equation}

Both approximations have truncation error ${\cal O}(h)$ linear in the step h.

More accurate is the central difference formula with error ${\cal O}(h^2)$

\begin{equation}
  f'(x) = \frac{f(x + h) - f(x - h)}{2 h}
\end{equation}

The purpose of this exercise is to calculate the numerical derivatives of a 
few simple functions as a function of the step $h$ and compare to 
the exact value calculated analytically. By going from very small to very large
steps $h$ (from $10^-{20}$ to $10^0$) one can also distinguish round off errors 
from truncation errors. In particular we will calculate the derivatives
of the following functions at the given point:

\begin{tabular}{c c c}
$x^2$ & at& $x=1$\\
$x^5$ & at &  $x=1$\\
$sin(x)$ & at  &$x=0$\\
$sin(x)$ & at  &$x=pi/4$\\
$sin(x)$ & at  &$x=pi/2$\\
\end{tabular}

Besides, the exercise is useful to learn other features of Matlab.
You want to write functions which calculate the forward, backward and central 
difference of a generic function. This means that the function itself has to be passed as argument. 
You have to do :
\begin{itemize}
\item create files fx2.m, fx5.m and fsin.m which calculate function and the 
exact derivatives for the functions $x^2$, $x^5$ and $sin(x)$.

\item create files forward\_difference.m ( etcetera) which calculate the forward
difference for a generic function fname. For instance:

{\tt function Forw\_Diff = Forward\_Diff(fname, x ,h)\\
\% Forward\_Diff approximates the derivative of a function f at point x\\
\% using the following two-points difference formula:\\
\% \\
\% derivative = ( f(x+h) - f(x) ) / h\\
\%\\
\% where h is the stepsize. The error is of order 1.\\
\\
Forw\_Diff = (feval(fname,x+h)-feval(fname,x))/h;}

In this example fname is a character string which gives the name of 
the function, {\tt feval} is a matlab 
function which evaluates the function fname in x. In the most recent
version of matlab one can refer to a function also with @fname where fname is 
the name of an m-file, e.g. feval(@fname,x). 

\item create a script derivatives.m which makes use of the previous files.
In this file you have to do several things

\begin{itemize}
\item read from the screen  the name of the 
functions of which you want to calculate the derivative and the value of the 
argument at which it has to be calculated. The name of the function is 
read as a string by using the command input with 's'. If 's' is not given 
matlab tries to evaluate the function and gives an error because the 
argument of the function is undefined. 
We also read the value x at which the derivative is evaluated.

{\tt  
fname = input('type in function of which derivative has to be \\
evaluated','s')\\
x = input ('type x, value at which the derivative will be evaluated')}
\item create a grid of values of h.
\item create vectors containing the values calculated with the
different approximations for all values of $h$ 

\item evaluate also the exact derivative for all values of $h$

\item calculate the relative error defined as
\begin{equation}
  \Delta = \frac {{f'}_{approx}(x) - f'(x)}{f'(x)}
\end{equation}

\item plot the calculated error with the different approximation, each 
with a different linetype.

\item functions that can be expressed as exponential or power law are 
conveniently plotted in logarithmic scale. Taking the logarithm of
the error $\Delta\sim C h^n$ where C is a constant and n is the order of the 
approximation, one obtains $ln(\Delta)=ln(C)+ nln(h)$ so that the slope 
in logarithmic scale gives directly the order of the approximation. 
Analogously one would obtain a straight line for an exponential function 
$y=b e^{ax}$ 
in a semilogarithmic plot since $ln(y)=ln(b)+ax$
\end{itemize}
\end{itemize}

\end{document}