\section{NUMERICAL DIFFERENTIATION}
The problem of computing derivatives is inherently unstable. To see
this, consider computing $f'(x_0)$ using the definition:
$$f'(x_0)=\lim_{h\rightarrow 0}{f(x_0+h)-f(x_0)\over h}.$$ In
floating-point arithmetic, as $h\to 0$, there will be serious
cancellation in the difference $f(x_0+h)-f(x_0)$, which is therefore
calculated with a large absolute error. Then the division by a small
$h$ results in an even larger absolute error. For $h$ small enough, the
computed number fl$\bigl[(f(x_0+h)-f(x_0))/h\bigr]$ is meaningless.
Intuitively, differentiation accentuates small errors, whereas
integration damps them out.

{\bf Ways to approximate} $f'(x)$:
\item{1)} Interpolate $f$ at $x_0$, $\ldots$, $x_n$ by a polynomial
$P_n$ and use $P_n'(x)$. The error in this case is
$$f^{(k)}(x)-P_n^{(k)}(x)={f^{(n-k+1)}(\zeta)\over (n-k+1)!}
\prod_{i=0}^{n-k} (x-\zeta_i),$$ where $x_i\le \zeta_i\le x_{i+k}$ and
$\zeta$ is in the smallest interval containing $x_0$, $\ldots$, $x_n$,
$x$.
\item{2)} Interpolate $f$ at $x_0<\cdots<x_n$ with a complete cubic
spline $S(x)$ and use $S'(x)$. Here the error satisfies $\left|
f^{(k)}(x) -S^{(k)}(x)\right| ={\cal O}(h^{4-k})$, $k=0$, 1, 2, 3,
$h=\max(x_{i+1}-x_i)$.
\item{3)} Use a {\it finite difference approximation}. Examples of
difference approximations to $f'(x)$ are:

Forward difference: $\displaystyle f'(x)={f(x+h)-f(x)\over h}+{\cal
O}(h)$.

Backward difference: $\displaystyle f'(x)={f(x)-f(x-h)\over h}+{\cal
O}(h)$.

Centered difference: $\displaystyle f'(x)={f(x+h)-f(x-h)\over h}+{\cal
O}(h^2)$.

Monday afternoon difference: $$f'(x)=
{-4f(x+2h)+18f(x+h)-21f(x)+10f(x-h)-3f(x-2h)\over 6h}+{\cal O}(h^3).$$

It is possible to get a difference approximation $\displaystyle\sum_{i=-k}^k
c_if(x+ih)$ to any derivative $f^{(m)}(x)$ with any accuracy ${\cal
O}(h^n)$, but generally high order formulas are not worthwhile. For
fixed $x$ and small $h>0$
let $\displaystyle D_hf={f(x+h)-f(x)\over h}$. Then
$$\displaylines{
\eqalign{f'(x)&=D_hf+a_1h+a_2h^2+\cdots+{\cal O}(h^k)\cr
f'(x)&=D_{qh}f+a_1(qh)+a_2(qh)^2+\cdots+{\cal O}(h^k)\cr}
\qquad\Longrightarrow\cr
f'(x)={qD_hf-D_{qh}f\over q-1}+a_2^*(qh)^2+a_3^*(qh)^3+\cdots+ {\cal
O}(h^k)=D_{qh}^*f+a_2^*(qh)^2+\cdots+{\cal O}(h^k).\cr}$$
Thus $D_{qh}^*$ is
a better approximation than $D_hf$ to $f'(x)$. $q$ can be anything, but
is usually taken to be 2. Now applying the same process to $D^*$ gives
$$\displaylines{
\eqalign{f'(x)&=D_{qh}^*f+a_2^*(qh)^2+a_3^*(qh)^3+\cdots\cr
f'(x)&=D_{q^2h}^*f+a_2^*(q^2h)^2+a_3^*(q^2h)^3+\cdots}\qquad
\Longrightarrow\cr
f'(x)={q^2D_{qh}^*f-D_{q^2h}^*f\over q^2-1}+a_3^{**}(q^2h)^3+\cdots
=D_{q^2h}^{**}f+{\cal O}(h^3).\cr}$$
This process can be continued if $f$ is
sufficiently differentiable to get better and better approximations (at
least higher order approximations). The resemblance of this process to
Romberg integration is clear, and these are examples of a general
procedure called {\it extrapolation to the limit}.

\noindent{\bf Theorem 1.} Let $A(y)=a_0+a_1y+a_2y^2+\cdots+a_ky^k+{\cal
O}(y^{k+1})$, $y>0$. Let $0<r<1$, and define a triangle of numbers
$A_{m,n}$, $n=0,1,\ldots,m$; $m=0,1,\ldots$ by
$$A_{m,0}=A(r^my),\qquad
A_{m,n}={A_{m,n-1}-r^nA_{m-1,n-1}\over 1-r^n}, \qquad n=1,\ldots,m.$$
Then $A_{m,n}=a_0+{\cal O}\left((r^my)^{n+1}\right)$, $n\le k$, and if
$a_n\ne0$, the $n$th column of the triangle converges to $a_0$ faster
than the $(n-1)$st column.

Proof. Henrici, {\it Elements of Numerical Analysis}, 1964.

\noindent{\bf Operator Calculus.} Often formulas can be derived by
purely {\it formal\/} manipulations, and then verified later by other
means. Consider the following linear operators on, say,
$C^\infty(-\infty,\infty)$: $$\eqalign{Ef(x)&=f(x+h)\cr \Delta 
f(x)&=f(x+h)-f(x)\cr Df(x)&=f'(x)\cr
\delta f(x)&=f\left(x+{h\over2}\right)-f\left(x-{h\over2}\right)\cr 
\mu f(x)&={1\over2}\left[f\left(x+{h
\over 2}\right)+f\left(x-{h\over2}\right)\right]\cr 
\nabla f(x)&=f(x)-f(x-h)\cr}$$
The Taylor expansion $\displaystyle 
f(x+h)=f(x)+(hD)f(x)+{(hD)^2\over 2!}f(x)+\cdots$
gives $$Ef(x)=e^{hD}f(x)\Rightarrow E=e^{hD}\Rightarrow
1+\Delta=e^{hD}\Rightarrow hD=\ln(1+\Delta)=\Delta-{1\over2}\Delta^2+{1
\over3}\Delta^3-\cdots.$$
This differentiation formula is valid if $f$ is
a polynomial or exponential. Another example is
$$\displaylines{\delta
f(x) =f\left(x+{1\over2}h\right)-f\left(x-{1\over2}h\right)
=\left(e^{hD/2}-e^{-hD/2}\right)f(x) =2\sinh\left({hD
\over 2}\right)f(x) \qquad
\Longrightarrow\cr \delta=2\sinh {hD\over2}.\cr}$$ 
Using $$ x=\sinh
x -{1\over2}{\sinh^3x\over3} +{1\over2}{3\over4}{\sinh^5x\over5} -{1\over
2}{3\over4}{5\over6}{\sinh^7x\over7} +\cdots$$
gives $$ hD =\delta-
{\delta^3\over24} +{3\delta^5\over640} -{5\delta^7\over7168} +\cdots=
\mu\delta -{\mu\delta^3\over6} +{\mu\delta^5\over30} -\cdots$$ since
$$1= \mu\left(1 +{1\over4}\delta^2\right)^{-1/2} =\mu\left(1
-{\delta^2\over8} +{3\delta^4
\over128} -{5\delta^6\over1024} +\cdots\right).$$

Here is a useful table of operator relationships:

\smallskip\centerline{\vbox{\tabskip=5pt\halign{&\hfil$#$\hfil\cr
\tablerule
\strut&E&\Delta&\delta&\nabla&D\cr
\noalign{\hrule\smallskip}
E&E&1+\Delta&1+{1\over2}\delta^2+\delta\sqrt{1+{1\over4}\delta^2}&
 {1\over1-\nabla}&e^{hD}\cr
\Delta&E-1&\Delta&\delta\sqrt{1+{1\over4}\delta^2}+{1\over2}\delta^2&
 {\nabla\over1-\nabla}&e^{hD}-1\cr
\delta&E^{1/2}-E^{-1/2}&\Delta(1+\Delta)^{-1/2}&\delta&
\nabla(1-\nabla)^{-1/2}&2\sinh{1\over2}hD\cr
\nabla&1-E^{-1}&{\Delta\over1+\Delta}&\delta\sqrt{1+{1\over4}
\delta^2}-{1\over2}\delta^2&\nabla&1-e^{-hD}\cr
D&{1\over h}\ln E&{1\over h}\ln(1+\Delta)&{2\over
h}\sinh^{-1}{1\over2}\delta&{1\over h}\ln\left({1\over1-\nabla}\right)&
D\cr
\mu&{1\over2}(E^{1/2}+E^{-1/2})&\left(1+{1\over2}\Delta\right)
(1+\Delta)^{-1/2}&\sqrt{1+{1\over4}\delta^2}&\left(1-{1\over2}\nabla
\right)(1-\nabla)^{-1/2}&\cosh{1\over2}hD\cr
\noalign{\smallskip\hrule}
}}}

\bigskip\centerline{DIFFERENCE EQUATIONS.}
A $k$th order difference equation has the form
$y_n=f(n,y_{n-1},y_{n-2}, \ldots,y_{n-k})$. A solution to a difference
equation is a sequence $\{y_n\}$, which is uniquely determined by
specifying $y_0,\ldots,y_{k-1}$. A $k$th order linear difference
equation with constant coefficients has the form $$ y_n +a_1y_{n-1}
+a_2y_{n-2} +\cdots +a_ky_{n-k} =b_n,\quad a_k\ne0.\leqno(1)$$ The
theory of difference equations closely parallels that of differential
equations. It is particularly useful for analyzing numerical methods.
The following theorems for 2nd order linear homogeneous difference
equations show the resemblance to differential equations. Consider the
difference equation $$ y_n +a_1y_{n-1} +a_2y_{n-2} =0,\quad
a_2\ne0.\leqno(2)$$ Let $Y^{(1)}=\{y_n^{(1)}\}$ and
$Y^{(2)}=\{y_n^{(2)}\}$ be solutions of (2). The {\it Wronskian\/} of
$Y^{(1)}$ and $Y^{(2)}$ at $n$ is $$w_n
=\left|\matrix{y_n^{(1)}&y_n^{(2)}\cr y_{n-1}^{(1)}&y_{n-1}^{(2)}\cr
}\right|.$$The characteristic equation of (2) is $t^2+a_1t+a_2=0$.

\noindent{\bf Theorem 2.} Let $Y^{(1)}$ and $Y^{(2)}$ be solutions of
(2). Then $c_1Y^{(1)}+c_2Y^{(2)}$ is also a solution.

\noindent{\bf Lemma 1.} A solution of (2) is uniquely determined by its
values at two consecutive integers.

\noindent{\bf Theorem 3.} Let $Y^{(1)}$ and $Y^{(2)}$ be solutions of
(2). Then every solution of (2) is a linear combination of
$Y^{(1)}$ and $Y^{(2)}\Leftrightarrow w_n\ne0$ $\forall n$.

\noindent{\bf Theorem 4.} Let $Y^{(1)}$ and $Y^{(2)}$ be solutions of
(2). Then $Y^{(1)}$ and $Y^{(2)}$ are independent $\Leftrightarrow
w_n\ne0$ for some $n$.

\noindent{\bf Theorem 5.} Let $Y^{(1)}$ and $Y^{(2)}$ be solutions of
(2). Then the Wronskian $w_n$ is either always zero or never zero.

\noindent{\bf Theorem 6.} Let $Y^{(1)}$ and $Y^{(2)}$ be solutions of
(2). Then every solution of (2) has the form $c_1Y^{(1)} +c_2Y^{(2)}
\Leftrightarrow w_n\ne0$ for some $n\Leftrightarrow Y^{(1)}$ and
$Y^{(2)}$ are independent.

\noindent{\bf Theorem 7.} Let $u_1, u_2$ be the roots of the
characteristic equation of (2). Then the general solution $\{y_n\}$ of
(2) is $$y_n=\cases{c_1u_1^n+c_1u_2^n,&if $u_1\ne u_2$,\cr
(c_1+c_2n)u_1^n, &if $u_1=u_2$.\cr}$$

All of the above theorems generalize to $k$th order difference equations
in the obvious way, e.g.,

\noindent{\bf Theorem 8.} Let $u_1,\ldots,u_m$ be the roots, with
multiplicity $r_1,\ldots,r_m$, of the characteristic equation of the
difference equation $y_n+a_1y_{n-1}+\cdots+a_ky_{n-k}=0$, $a_k\ne 0$.
Then the general solution $\{y_n\}$ is given by $$y_n =P_1(n)u_1^n
+P_2(n)u_2^n+ \cdots+ P_m(n)u_m^n,$$ where $P_i(n)$ is an arbitrary
polynomial of degree $r_i-1$ for $i=1$, $\ldots$, $m$.

A $k$th order difference equation can be reduced to a system of first
order difference equations. The nonhomogeneous first order linear system
has the form $Y_{n+1}=B_nY_n+X_n$, where $\{X_n\}$ is a given sequence
of vectors. The solution is $$Y_n=P_{n0}Y_0+\sum_{j=1}^nP_{nj}X_{j-1},$$
where $P_{nj}=B_{n-1}B_{n-2}\cdots B_j$, $P_{nn}=I$.