%
%
%

\section{Numerical Differentiation}   \label{sec_int01}
\begin{frame}
\frametitle{Numerical Differentiation}

\begin{eqnarray*}
    f^\prime(x_0) = \lim_{h \to 0} \frac{f(x_0+h) - f(x_0)}{h}.
\end{eqnarray*}

\begin{block}{Question}
    How accurate is
    \begin{eqnarray*}
        \frac{f(x_0+h) - f(x_0)}{h}?
    \end{eqnarray*}
\end{block}



Suppose a given function $f$ has continuous first derivative and
$f''$ exists. From Taylor's theorem
%
\[
   f(x + h) = f(x) + f'(x) h + \frac{1}{2} f''(\xi) h^{2},
\]
%
where $\xi$ is between $x$ and $x+h$, one has
%
\[
   f'(x) = \frac{f(x+h) - f(x)}{h} - \frac{h}{2} f''(\xi)
   = \frac{f(x+h) - f(x)}{h} + O(h).
\]
\end{frame}
%

\begin{frame}
Hence it is reasonable to use the approximation
%
\begin{eqnarray*}
   f'(x) \approx
   \frac{f(x+h) - f(x)}{h}
\end{eqnarray*}
which is called forward finite difference, and the error involved is
%
\begin{eqnarray*}
   |e| = \frac{h}{2} |f''(\xi)| \leq \frac{h}{2}
   \max_{t \in (x,x+h)} |f''(t)|.
\end{eqnarray*}


Similarly one can derive the backward finite difference
approximation
%
\begin{equation}
   f'(x) \approx \frac{f(x) - f(x-h)}{h}
\end{equation}
%
which has the same order of truncation error as the forward
finite difference scheme.
\end{frame}

\begin{frame}

The forward difference is an $O(h)$ scheme. An $O(h^{2})$ scheme
can also be derived from the Taylor's theorem
%
\begin{eqnarray*}
   f(x + h)
   & = &
   f(x) + f'(x) h + \frac{1}{2} f''(x) h^{2}
   + \frac{1}{6} f'''(\xi_{1}) h^{3}   \\
   f(x-h)
   & = &
   f(x) - f'(x) h + \frac{1}{2} f''(x) h^{2}
   - \frac{1}{6} f'''(\xi_{2}) h^{3},
\end{eqnarray*}
%
where $\xi_{1}$ is between $x$ and $x+h$ and $\xi_{2}$ is between
$x$ and $x-h$. Hence
\[
   f(x+h) - f(x-h) = 2 f'(x) h +
   \frac{1}{6} [f'''(\xi_{1}) + f'''(\xi_{2})] h^{3}
\]
%
and
%
\[
   f'(x) = \frac{f(x+h) - f(x-h)}{2h}
   - \frac{1}{12} [f'''(\xi_{1}) + f'''(\xi_{2})] h^{2}
\]
%
Let
\[
    M = \max_{z \in [x-h,x+h]} f'''(z)
    \quad \text{and} \quad
    m = \min_{z \in [x-h,x+h]} f'''(z).
\]
\end{frame}

%
\begin{frame}
If $f^{\prime\prime\prime}$ is continuous on $[x-h,x+h]$, then by
the intermediate value theorem, there exists $\xi \in [x-h,x+h]$
such that
\begin{eqnarray*}
    f^{\prime\prime\prime}(\xi) = \frac{1}{2}
    [f^{\prime\prime\prime}(\xi_{1}) +
    f^{\prime\prime\prime}(\xi_{2})].
\end{eqnarray*}
Hence
%
\begin{eqnarray*}
   f^{\prime}(x) = \frac{f(x+h) - f(x-h)}{2h} - \frac{1}{6}
   f^{\prime\prime\prime}(\xi) h^{2}
   =  \frac{f(x+h) - f(x-h)}{2h} +  O(h^{2}).
\end{eqnarray*}
%
This is called center difference approximation and the truncation
error is
%
\[
   |e| = \frac{h^{2}}{6} f^{\prime\prime\prime}(\xi)
\]



Similarly, we can derive an $O(h^{2})$ scheme from Taylor's theorem
for $f^{\prime\prime}(x)$
\begin{eqnarray*}
   f^{\prime\prime}(x) = \frac{f(x+h) - 2 f(x) + f(x-h)}{h^{2}}
   - \frac{1}{12} f^{(4)}(\xi) h^{2},
\end{eqnarray*}
where $\xi$ is between $x-h$ and $x+h$.

\end{frame}

\begin{frame}
\frametitle{Polynomial Interpolation Method}


%A general approach for numerical differentiation is to use
%polynomial interpolation.


Suppose that $(x_{0},f(x_{0}))$, $(x_{1},f(x_{1}))$ $\cdots, (x_{n},
f(x_{n}))$ have been given, we apply the Lagrange polynomial
interpolation scheme to derive
%
\begin{eqnarray*}
   P(x) = \sum_{i=0}^{n} f(x_{i}) L_{i}(x),
\end{eqnarray*}
%
where
%
\begin{eqnarray*}
   L_{i} (x) = \prod_{j=0,j \neq i}^{n}
   \frac{x - x_{j}}{x_{i} - x_{j}}.
\end{eqnarray*}
%
Since $f(x)$ can be written as
%
\begin{eqnarray*}
   f(x) = \sum_{i=0}^{n} f(x_{i}) L(x)
   + \frac{1}{(n+1)!} f^{(n+1)}(\xi_{x}) w(x),
\end{eqnarray*}
%
where
%
\begin{eqnarray*}
   w(x) = \prod_{j=0}^{n} (x - x_{j}),
\end{eqnarray*}
\end{frame}
%

\begin{frame}
we have,
%
\begin{eqnarray*}
   f'(x) &=& \sum_{i=0}{n} f(x_{i}) L_{i}'(x) +
   \frac{1}{(n+1)!} f^{(n+1)}(\xi_{x}) w'(x) \\
   &+& \frac{1}{(n+1)!} w(x) \frac{d}{dx} f^{(n+1)} (\xi_{x}).
\end{eqnarray*}
%
Note that
%
\begin{eqnarray*}
   w'(x) = \sum_{j=0}^{n} \prod_{i=0,i \neq j}^{n} (x - x_{i}).
\end{eqnarray*}
%
Hence a reasonable approximation for the first derivative of
$f$ is
%
\begin{eqnarray*}
   f'(x) \approx \sum_{i=0}^{n} f(x_{i}) L_{i}'(x).
\end{eqnarray*}
%
When $x = x_{k}$ for some $0 \leq k \leq n$,
%
\[
    w(x_{k}) = 0 \quad \text{and} \quad
    w'(x_{k}) = \prod_{i=0, i \neq k}^{n} (x_{k} - x_{i}).
\]
\end{frame}

%
\begin{frame}
Hence
%
\begin{equation}
   f'(x_{k}) = \sum_{i=0}^{n} f(x_{i}) L_{i}' (x_{k})
   + \frac{1}{(n+1)!} f^{(n+1)}(\xi_{x})
   \prod_{i=0, i \neq k}^{n} (x_{k} - x_{i}),
\end{equation}
which is called an $(n+1)$-point formula to approximate
$f^\prime(x)$.

$\bullet$ Three Point Formulas

Since
\begin{eqnarray*}
    L_0(x) = \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}
\end{eqnarray*}
we have
\begin{eqnarray*}
   L_0^{\prime}(x) = \frac{2x - x_1 - x_2}{(x_0-x_1)(x_0-x_2)}.
\end{eqnarray*}
Similarly,
\begin{eqnarray*}
   L_1^{\prime}(x) = \frac{2x - x_0 - x_2}{(x_1-x_0)(x_1-x_2)} \
   \mbox{ and } \    L_2^{\prime}(x) = \frac{2x - x_0 -
   x_1}{(x_2-x_0)(x_2-x_1)}.
\end{eqnarray*}
\end{frame}

\begin{frame}
Hence
\begin{eqnarray*}
    f^\prime(x_j) &=& f(x_0) \left[ \frac{2x_j - x_1 -
    x_2}{(x_0-x_1)(x_0-x_2)} \right] + f(x_1) \left[ \frac{2x_j - x_0 -
    x_2}{(x_1-x_0)(x_1-x_2)} \right] \\
    &+& f(x_2) \left[ \frac{2x_j - x_0 -
   x_1}{(x_2-x_0)(x_2-x_1)} \right] + \frac{1}{6} f^{(3)}(\xi_j)
   \prod_{k= 0, k \neq j}^2 ( x_j - x_k),
\end{eqnarray*}
for each $j = 0, 1, 2$. Assume that
\begin{eqnarray*}
   x_1 = x_0 + h \mbox{ and } \ x_2 = x_0 + 2 h, \ \mbox{ for some }
   \ h \neq 0.
\end{eqnarray*}
Then
\begin{eqnarray*}
    f^{\prime}(x_0) &=& \frac{1}{h} \left[ -\frac{3}{2} f(x_0) + 2
    f(x_1) - \frac{1}{2} f(x_2) \right] + \frac{h^2}{3}
    f^{(3)}(\xi_0), \\
    f^{\prime}(x_1) &=& \frac{1}{h} \left[ - \frac{1}{2} f(x_0) +
    \frac{1}{2} f(x_2) \right] - \frac{h^2}{6} f^{(3)}(\xi_1), \\
    f^{\prime}(x_2) &=& \frac{1}{h} \left[ \frac{1}{2} f(x_0) - 2
    f(x_1) + \frac{3}{2} f(x_2) \right] + \frac{h^2}{3}
    f^{(3)}(\xi_2).
\end{eqnarray*}
\end{frame}

\begin{frame}
That is
\begin{eqnarray}
    f^{\prime}(x_0) &=& \frac{1}{h} \left[ -\frac{3}{2} f(x_0) + 2
    f(x_0+h) - \frac{1}{2} f(x_0+2h) \right] + \frac{h^2}{3}
    f^{(3)}(\xi_0), \nonumber \\
    f^{\prime}(x_0+h) &=& \frac{1}{h} \left[ - \frac{1}{2} f(x_0) +
    \frac{1}{2} f(x_0+2h) \right] - \frac{h^2}{6} f^{(3)}(\xi_1), \label{eq:4.4_01} \\
    f^{\prime}(x_0+2h) &=& \frac{1}{h} \left[ \frac{1}{2} f(x_0) - 2
    f(x_0+h) + \frac{3}{2} f(x_0+2h) \right] + \frac{h^2}{3}
    f^{(3)}(\xi_2). \label{eq:4.4_02}
\end{eqnarray}
Using the variable substitution $x_0$ for $x_0+h$ and $x_0+2h$ in
\eqref{eq:4.4_01} and \eqref{eq:4.4_02}, respectively, we have
\begin{eqnarray}
    f^{\prime}(x_0) &=& \frac{1}{2h} \left[ -3 f(x_0) + 4
    f(x_0+h) - f(x_0+2h) \right] + \frac{h^2}{3}
    f^{(3)}(\xi_0), \label{eq:4.4_03} \\
    f^{\prime}(x_0) &=& \frac{1}{2h} \left[ - f(x_0-h) +
    f(x_0+h) \right] - \frac{h^2}{6} f^{(3)}(\xi_1),  \nonumber \\
    f^{\prime}(x_0) &=& \frac{1}{2h} \left[ f(x_0-2h) - 4
    f(x_0-h) + 3 f(x_0) \right] + \frac{h^2}{3}
    f^{(3)}(\xi_2). \label{eq:4.4_04}
\end{eqnarray}
Note that \eqref{eq:4.4_04} can be obtained from \eqref{eq:4.4_03}
by replacing $h$ with $-h$.
\end{frame}

\begin{frame}
$\bullet$ Five-point Formulas

\begin{eqnarray*}
    f^{\prime}(x_0) &=& \frac{1}{12h} \left[ f(x_0-2h) - 8 f(x_0-h) +
    8 f(x_0+h) - f(x_0+2h) \right] \\
    &+& \frac{h^4}{30} f^{(5)}(\xi),
\end{eqnarray*}
where $\xi \in (x_0 - 2h, x_0+2h)$ and
\begin{eqnarray*}
    f^{\prime}(x_0) &=& \frac{1}{12h} \left[ -25 f(x_0) + 48 f(x_0+h)
    - 36 f(x_0+2h) \right. \\
    &+& \left. 16 f(x_0+3h) - 3 f(x_0+4h) \right] + \frac{h^4}{5}
    f^{(5)}(\xi),
\end{eqnarray*}
where $\xi \in (x_0, x_0+4h)$.
\end{frame}

\begin{frame}
\frametitle{Round-off Error}

Consider
\begin{eqnarray*}
    f^{\prime}(x_0) &=& \frac{1}{2h} \left[ - f(x_0-h) +
    f(x_0+h) \right] - \frac{h^2}{6} f^{(3)}(\xi_1),
\end{eqnarray*}
where $\frac{h^2}{6} f^{(3)}(\xi_1)$ is called truncation error. Let
$\tilde{f}(x_0+h)$ and $\tilde{f}(x_0-h)$ be the computed values of
$f(x_0+h)$ and $f(x_0-h)$, respectively. Then
\begin{eqnarray*}
    f(x_0+h) = \tilde{f}(x_0+h) + e(x_0+h)
\end{eqnarray*}
and
\begin{eqnarray*}
    f(x_0-h) = \tilde{f}(x_0-h) + e(x_0-h).
\end{eqnarray*}
Therefore, the total error in the approximation
\begin{eqnarray*}
    f^{\prime}(x_0) - \frac{\tilde{f}(x_0+h) - \tilde{f}(x_0-h)}{2h}
    = \frac{e(x_0+h) - e(x_0-h)}{2h} - \frac{h^2}{6} f^{(3)}(\xi_1)
\end{eqnarray*}
is due in part to round-off error and in the part to truncation
error.
\end{frame}

\begin{frame}
Assume that
\begin{eqnarray*}
    | e(x_0 \pm h) | \leq \varepsilon \ \mbox{ and } \ | f^{(3)}(
    \xi_1 ) | \leq M.
\end{eqnarray*}
Then
\begin{eqnarray*}
    \left| f^{\prime}(x_0) - \frac{\tilde{f}(x_0+h) -
    \tilde{f}(x_0-h)}{2h} \right|
    \leq \frac{\varepsilon}{2h} + \frac{h^2}{6} M \equiv e(h).
\end{eqnarray*}
Note that $e(h)$ attains its minimum at $h =
\sqrt[3]{3\varepsilon/M}$.

\end{frame}

\section{Richardson Extrapolation Method}
\begin{frame}
\frametitle{Richardson's Extrapolation}

Suppose $\forall h \neq 0$ we have a formula $N_1(h)$ that
approximates an unknown value $M$
\begin{eqnarray}
     M - N_1(h) = K_1 h + K_2 h^2 + K_3 h^3 + \cdots, \label{eq:Ric_1}
\end{eqnarray}
for some unknown constants $K_1, K_2, K_3, \ldots$. If $K_1 \neq 0$,
then the truncation error is $O(h)$. For example,
\begin{eqnarray*}
   f^\prime(x) - \frac{f(x+h)-f(x)}{h}  =
   - \frac{f''(x)}{2!} h
   - \frac{f'''(x)}{3!}h^{2}  - \frac{f^{(4)}(x)}{4!}h^{3}
   -  \cdots.
\end{eqnarray*}

\begin{block}{Goal}
    Find an easy way to produce formulas with a higher-order
    truncation error.
\end{block}

Replacing $h$ in \eqref{eq:Ric_1} by $h/2$, we have
\begin{eqnarray}
    M = N_1 \left( \frac{h}{2} \right) + K_1 \frac{h}{2} + K_2
    \frac{h^2}{4} + K_3 \frac{h^3}{8} + \cdots. \label{eq:Ric_2}
\end{eqnarray}

\end{frame}

\begin{frame}
Subtracting \eqref{eq:Ric_1} with twice \eqref{eq:Ric_2}, we get
\begin{eqnarray}
    M = N_2(h) - \frac{K_2}{2} h^2 - \frac{3K_3}{4} h^3 - \cdots,
    \label{eq:Ric_3}
\end{eqnarray}
where
\begin{eqnarray*}
    N_2(h) = 2 N_1 \left( \frac{h}{2} \right) - N_1(h) = N_1 \left( \frac{h}{2}
    \right) + \left[ N_1 \left( \frac{h}{2} \right) - N_1(h)
    \right],
\end{eqnarray*}
which is an $O(h^2)$ approximation formula.

Replacing $h$ in \eqref{eq:Ric_3} by $h/2$, we get
\begin{eqnarray}
    M = N_2\left( \frac{h}{2} \right) - \frac{K_2}{8} h^2 -
    \frac{3K_3}{32} h^3 - \cdots. \label{eq:Ric_4}
\end{eqnarray}
Subtracting \eqref{eq:Ric_3} from 4 times \eqref{eq:Ric_4} gives
\begin{eqnarray*}
    3M = 4N_2\left( \frac{h}{2} \right) - N_2(h) + \frac{3K_3}{8} h^3
    + \cdots,
\end{eqnarray*}
which implies that
\begin{eqnarray*}
    M = \left[ N_2 \left( \frac{h}{2} \right) + \frac{N_2(h/2) -
    N_2(h)}{3} \right] + \frac{K_3}{8} h^3 + \cdots \equiv N_3(h) +
    \frac{K_3}{8} h^3 + \cdots.
\end{eqnarray*}
\end{frame}

\begin{frame}
Using induction, $M$ can be approximated by
\begin{eqnarray*}
    M = N_m(h) + O(h^m),
\end{eqnarray*}
where
\begin{eqnarray*}
    N_m(h) = N_{m-1}\left( \frac{h}{2} \right) + \frac{N_{m-1}(h/2)
    - N_{m-1}(h)}{2^{m-1} - 1}.
\end{eqnarray*}

{\bf Centered difference formula}. From the Taylor's theorem
%
\begin{eqnarray*}
   f(x+h) & = &
   f(x) + h f'(x) + \frac{h^{2}}{2!} f''(x)
   + \frac{h^{3}}{3!} f'''(x) + \frac{h^{4}}{4!} f^{(4)}(x)
   + \frac{h^{5}}{5!} f^{(5)}(x) + \frac{h^{6}}{6!} f^{(6)}(x)
   + \cdots  \\
   f(x-h) & = &
   f(x) - h f'(x) + \frac{h^{2}}{2!} f''(x)
   - \frac{h^{3}}{3!} f'''(x) + \frac{h^{4}}{4!} f^{(4)}(x)
   - \frac{h^{5}}{5!} f^{(5)}(x) + \frac{h^{6}}{6!} f^{(6)}(x)
   + \cdots
\end{eqnarray*}
%
we have
%
\[
   f(x+h) - f(x-h) = 2 h f'(x) + \frac{2 h^{3}}{3!} f'''(x)
   + \frac{2 h^{5}}{5!} f^{(5)}(x) + \cdots,
\]
\end{frame}

%
\begin{frame}
and, consequently,
%
\begin{eqnarray}
   f'(x_0) &=& \frac{f(x_0+h) - f(x_0-h)}{2h} - \left[
   \frac{h^{2}}{3!} f'''(x_0) + \frac{h^{4}}{5!} f^{(5)}(x_0)
   + \cdots \right], \nonumber \\
   &\equiv& N_1(h)- \left[
   \frac{h^{2}}{3!} f'''(x_0) + \frac{h^{4}}{5!} f^{(5)}(x_0)
   + \cdots \right]. \label{eq:Ric_5}
\end{eqnarray}
Replacing $h$ in \eqref{eq:Ric_5} by $h/2$ gives
\begin{eqnarray}
    f^\prime(x_0) = N_1\left( \frac{h}{2} \right) - \frac{h^2}{24}
    f^{\prime\prime\prime}(x_0) - \frac{h^4}{1920} f^{(5)}(x_0) -
    \cdots. \label{eq:Ric_6}
\end{eqnarray}
Subtracting \eqref{eq:Ric_5} from 4 times \eqref{eq:Ric_6} gives
\begin{eqnarray*}
    f^\prime(x_0) = N_2(h) + \frac{h^4}{480} f^{(5)}(x_0) + \cdots,
\end{eqnarray*}
where
\begin{eqnarray*}
    N_2(h) = \frac{1}{3} \left[ 4 N_1\left( \frac{h}{2} \right) -
    N_1(h) \right] = N_1\left( \frac{h}{2} \right) + \frac{N_1(h/2)
    - N_1(h)}{3}.
\end{eqnarray*}

\end{frame}

\begin{frame}
%
In general,
\begin{eqnarray*}
    f^\prime(x_0) = N_j(h) + O(h^{2j})
\end{eqnarray*}
with
\begin{eqnarray*}
    N_j(h) = N_{j-1} \left( \frac{h}{2} \right) + \frac{N_{j-1}(h/2)
    - N_{j-1}(h)}{4^{j-1} - 1}.
\end{eqnarray*}

\begin{example}
    Suppose that $x_0 = 2.0$, $h = 0.2$ and $f(x) = x e^x$. Compute
    an approximated value of $f^\prime(2.0) = 22.16716829679195$ to six decimal
    places.
\end{example}
{\bf Solution}. By centered difference formula, we have
\begin{eqnarray*}
    N_1(0.2) &=& \frac{f(2.0+0.2) - f(2.0-0.2)}{2h} = 22.414160,  \\
    N_1(0.1) &=& \frac{f(2.0+0.1) - f(2.0-0.1)}{h} = 22.228786.
\end{eqnarray*}
\end{frame}

\begin{frame}
It implies that
\begin{eqnarray*}
    N_2(0.2) = N_1(0.1) + \frac{N_1(0.1) - N_1(0.2)}{3} = 22.166995
\end{eqnarray*}
which does not have six decimal digits. Adding $N_1(0.05) =
22.182564$, we get
\begin{eqnarray*}
    N_2(0.1) = N_1(0.05) + \frac{N_1(0.05) - N_1(0.1)}{3} =
    22.167157
\end{eqnarray*}
and
\begin{eqnarray*}
    N_3(0.2) = N_2(0.1) + \frac{N_2(0.1) - N_2(0.2)}{15} = 22.167168
\end{eqnarray*}
which contains six decimal digits.
\end{frame}

\begin{frame}
\begin{tabular}[c]{llll} \hline
   \multicolumn{1}{c}{{\textcolor{blue}{$O(h)$}}} &
      \multicolumn{1}{c}{{\textcolor{blue}{$O(h^2)$}}} &
      \multicolumn{1}{c}{{\textcolor{blue}{$O(h^3)$}}} &
      \multicolumn{1}{c}{{\textcolor{blue}{$O(h^4)$}}} \\ \hline
   {\textcolor{red}{1}}: $N_1(h) = N(h)$ \\
   {\textcolor{red}{2}}: $N_1(h/2) = N(h/2)$ & {\textcolor{red}{3}}: $N_2(h)$ \\
   {\textcolor{red}{4}}: $N_1(h/4) = N(h/4)$ & {\textcolor{red}{5}}: $N_2(h/2)$ &
      {\textcolor{red}{6}}: $N_3(h)$ \\
   {\textcolor{red}{7}}: $N_1(h/8) = N(h/8)$ & {\textcolor{red}{8}}: $N_2(h/4)$ &
      {\textcolor{red}{9}}: $N_3(h/2)$ & {\textcolor{red}{10}}: $N_4(h)$
   \\ \hline
\end{tabular}
\end{frame}
%
