NUMERICAL DIFFERENTIATION

Recall that

$$f'(x_0)=\lim_{h\to 0}\frac{f(x_0+h)-f(x_0)}{h}.$$ (39)

Assume

$$\begin{eqnarray*} x_0\in (a,b),\quad f\in C^2[a,b]\\ x_1=x_0+h, \quad \mbox {with } h \ne 0\quad \mbox {small so that }x_1\in [a,b]. \end{eqnarray*}$$

$$f(x_0+h)\approx f(x_0)+hf'(x_0). \quad \mbox { More precisely}$$

$$f(x_0+h)=f(x_0)+hf'(x_0)+\frac12 h^2f''(\xi)\mbox { with }\quad \xi\in (x_0, x_0+h)$$

$$\therefore f'(x_0)=\frac{f(x_0+h)-f(x_0)}{h}-\frac12 hf''(\xi)$$ (40)

Comparison of (40) and (41) shows that replacing the $$\lim_{h\to 0}$$ by a difference as $h\to 0$ introduces an error

$$\frac12 hf''(\xi)=O(h)\quad \equiv \mbox {{\em truncation error}}.$$

(41) is the forward difference when $h>0$ and backward difference when $h<0$.

Example $f(x)= \ln x\quad x_0=1.8$

$$\frac{f(1.8+h)-f(1.8)}{h}\quad h>0$$

used to find $f'(1.8)$ with error

$$\frac{\vert hf''(\xi)\vert}{2}=\frac{\vert h\vert}{2\xi^2}\le \frac{\vert h\vert}{2(1.8)^2} \mbox { where } 1.8<\xi< 1.8+h$$

$\nwarrow$ from $f''$

$h$	$f'$	$\frac{\vert h\vert}{2\xi^2}$
$0.1$	$0.5406722$	$0.0154321$
$0.01$	$0.5546180$	$0.0015432$
$0.001$	$0.55554013$	$0.0001543$

$\Box$

More general derivative formulas:

Suppose $\{x_0, x_1,\cdots, x_n\}\mbox{ are } n+1$ distinct numbers, $I$ is the interval and $f\in C^{n+1}(I)$. Using Lagrange polynomials $\ell$:

Express

$$f(x)=\sum^{n}_{k=0}f(x_k)\ell_k(x)+\frac{(x-x_0)\cdots (x-x_n)} {(n+1)!}f^{(n+1)}(\xi)$$

for $\xi(x)\in I.$

Differentiating:

$$\begin{eqnarray*} f'(x)=\sum^{n}_{k=0}f(x_k)\ell'_k(x)+\frac{d}{d x}\left[\frac... ...0-x_0)\cdots(x-x_n) \frac{d}{d x} \left(f^{(n+1)}(\xi(x))\right) \end{eqnarray*}$$

What's the truncation error? If $x=x_k$ then term involving $D_x[f^{n+1}(\xi)]=0$. Therefore,

$$\begin{eqnarray*} f'(x_k)=\sum^n_{j=0}f(x_j)\ell'_j(x_k)+ \frac{f^{(n+1)}(\xi(x_k))}{(n+1)!}\Pi^n_{j=0\atop j\ne k} (x_k-x_j) \end{eqnarray*}$$

which is equivalent to the $n+1$-point formula, since

$$\frac{\partial }{\partial x}\prod^n_{j=0}(x-x_j)=\sum^n_{j=... ...x { if } x=x_k \Rightarrow\prod^n_{j=0\atop j\ne k}(x_k-x_j).$$

In general, using more points leads to greater accuracy.

Price payed: more functional evaluations and possible round-off error. Thus, we must balance truncation error versus round-off error and speed benefits.

Most common: 3 and 5 point formulas.

$$\begin{eqnarray*} &&\ell_0(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} \quad \... ...2)}\\ \\ &&\ell'_2(x)=\frac{(2x-x_0-x_1)}{(x_2-x_0)(x_2-x_1)} \end{eqnarray*}$$

Therefore,

$$f'(x_j)=f(x_0)\left[\frac{2xj-x_1-x_2}{(x_0-x_1)(x_0-x_2)}\right]$$

$$+f(x_1)\left[\frac{2xj-x_0-x_2}{(x_1-x_0)(x_1-x_2)}\right]$$ (41)

$$+f(x_2)\left[\frac{2xj-x_0-x_1}{(x_2-x_0)(x_2-x_1)}\right]+\frac16 f^{(3)}(\xi_j) \frac{d}{dx}\prod^2_{i=0\atop i\ne j}(x_j-x_i)$$

for each $j=0, 1, 2$ and $\xi_j=\xi(x_j)$

Look at (42) on an evenly spaced grid:

$x_1=x_0+h$    $x_2=x_0+2h$,     $h\ne 0$

take $x_j=x_0\quad, \quad x_1=x_0+h,\quad \quad x_2=x_0+2L$

$$f'(x_0)=\frac{1}{h}\left[-\frac32 f(x_0)+2f(x_1)-\frac12 f(x_2)\right]+\frac{h^2}{3}f^{(3)}(\xi_0)$$

when $x_j=x_i$

$$f'(x_1)=\frac{1}{h}\left[-\frac12 f(x_0)+\frac12 f(x_2)\right] -\frac{h^2}{3}f^{(3)}(\xi_1)$$

when $x_j=x_2$

$$f'(x_2)=\frac{1}{h}\left[\frac12f(x_0)-2f(x_1)+\frac32 f(x_2)\right]+\frac{h^2}{3}f^{(3)}(\xi_2)$$

Now, since $x_1=x_0+h, \quad x_2=x_0+2h\Rightarrow$

$$\begin{eqnarray*} f'(x_0)=\frac{1}{h}\left[-\frac32 f(x_0)+2f(x_0+h) -\frac12f(... ...)-2f(x_0+h)+\frac32f(x_0+2h) \right]+\frac{h^2}{3}f^{(3)}(\xi_2) \end{eqnarray*}$$

We can translate:

$$\begin{eqnarray*} f'(x_0)=\frac{1}{2h}\left[-3f (x_0)+4f(x_0+h)-f(x_0+2h)\right... ...f(x_0-2h)-4f(x_0-h)+3f(x_0)\right] +\frac{h^2}{3}f^{(3)}(\xi_2) \end{eqnarray*}$$

$\therefore$ 2 formulas, really, since the last and first are equivalent by letting $h\rightarrow -h$.

3-point formulas : Recall that

$$f'(x_0)=\frac{1}{2h}[-3f(x_0)+4f(x_0+h)-f(x_0+2h)]+ \frac{h^2}{3}f^{(3)}(\xi_0)$$ (42)

$$f'(x_0)=\frac{1}{2h}[f(x_0+h)-f(x_0-h)]-\frac{h^2}{6}f^{(3)} (\xi_1)$$ (43)

where

$$\begin{eqnarray*} \xi_0\,\mbox {lies between }x_0 \mbox { and } x_0+2h\\ \xi_1 \, \mbox {lies between } (x_0-h) \mbox { and } (x_0+h) \end{eqnarray*}$$

(44) error $\sim \displaystyle \frac12$ (43) error, since data is examined on both sides of $x_0$.

There are 5-point formulas

$$f'(x_0)=\frac{1}{12h}[f(x_0-2h)-8f(x_0-h)+8f(x_0+h)-f(x_0+2h)]$$

$$+\frac{h^4}{30}f^{5}(\xi)$$ (44)

$$f'(x)=\frac{1}{12h}[-25f(x_0)+48f(x_0+h)-36f(x_0+2h)+16f(x_0+3h)$$

$$-3f(x_0+4h)]+\frac{h^4}{5}f^{5}(\xi)$$ (45)

(Can use $h<0$ or $h>0$ for left and right)     $x_0\le \xi\le x_0+4L$

(43) is useful at end of interval.

(46) is also useful formula at end of interval.

Higher-order derivatives?

Can use same procedure as above. Another way:

Example Take $f''$ Second derivative of $f$, for example, found by expanding

$$f(x_0+h)=f(x_0)+f'(x_0)h+f^{('')}(x_0)\frac{h^2}{2}+f^{(3)}(x_0) \displaystyle \frac{h^3}{6}+f^{''''}(\xi_{+h})\frac{h^4}{24}$$

$f(x_0-h)=f(x_0)-f'(x_0)h+f^{''}(x_0)\displaystyle \frac{h^2}{2}-f^{(3)}(x_0) \displaystyle \frac{h^3}{6}+f^{(4)}(\xi_{-h})\frac{h^4}{24}$

and add and solve for $f^{''}$:

$$f^{''}(x_0)=\frac{1}{h^2}\Big[f(x_0-h)-2f(x_0)+f(x_0+h)\Big] -\frac{h^2}{24}\Big[f^{(4)}(\xi_1)+f^{(4)}(\xi_{-1})\Big]$$

If $f^{(4)}$ is continued on $[x_0-h, x_0+h]$ the intermediate value theorem permits

$$\begin{eqnarray*} f''(x_0)=\frac{1}{h^2}[f(x_0+h)-2f(x_0)+f(x_0-h)] -\frac{h^2}{12}f^{(4)}(\xi)\\ \mbox {for some } x_0-h<\xi<x_0+h. \end{eqnarray*}$$

$\Box$

Suppose we take into account truncation and round off errors?

Take $f'$ case: Let $\tilde f$ be computed approximation to $f$

$$f'(x_0)=\frac{1}{2h}[f(x_0+h)-f(x_0-h)]-\frac{h^2}{6}f^{(3)} (\xi_1)$$

(i) in evaluating                     $\nearrow$                     $\nearrow$     we encounter round-off

$e(x_0+h)$ and $e(x_0-h)$

$\Rightarrow f(x_0\pm h)=\tilde f(x_0\pm h)+e(x_0\pm h)$

$$\begin{eqnarray*} % latex2html id marker 3309&\therefore f'(x_0)-\frac{\tilde... ...ad \mbox { and } \big\vert f^{(3)}(\xi_1)\big\vert<M\quad M>0\\ \end{eqnarray*}$$

$$\begin{eqnarray*} % latex2html id marker 3321&&\therefore \Big\vert f'(x_0)-\... ...ll can make error be dominated by round-off.} (See \ref{rndoff}) \end{eqnarray*}$$

In Figure 43 we show the error as a function of $h$ for the central difference approximation of a smooth function.

Figure 43: Error in the central difference approximation to the derivative $f'(x)$ as a function of $h$. The dots correspond to the error solely due to truncation error. Stars correspond to the the error from both truncation and roundoff.

It is clear from the star curve, which corresponds to both the truncation and roundoff error contributions, that for small $h$ the error is dominated by the roundoff and for larger $h$ it is dominated by the truncation error. The optimal value of $h$, which minimizes the error can be found by minimizing the expression for the error.

The table below shows the errors for different values of $h$. Note that for larger values of $h$ ($h> 10^{-6}$) the errors seem to be decreasing by a factor of approximately 100. This is due to the truncation error of the central difference method. As the method is of order $O(h^2)$ when we decrease $h$ by ten we decrease the error by 100. But as $h$ becomes smaller the round off error of the computer becomes more important and the error starts to increase with decreasing $\epsilon$. This increase is approximately by a factor of 10 with every decrease of $h$ by a factor of 10. This is because the $\frac{\epsilon}{h}$ term is dominating for small values of $h$. The value of $M$ given in this example is $f'''(x_0)=4804.6$; the value of $\epsilon$ is $2 \times f(x_0) \times \mbox{machine precision} = 2.4246 \times 10^{-14}$.

$h$	$\vert$Error$\vert$	Estimated Error
$10^{-1}$	$8.124453$	$8.007729$
$10^{-2}$	$8.008886 \times 10^{-2}$	$8.007729\times 10^{-2}$
$10^{-3}$	$8.007740 \times 10^{-4}$	$8.007729\times 10^{-4}$
$10^{-4}$	$8.007935 \times 10^{-6}$	$8.007971\times 10^{-6}$
$10^{-5}$	$8.172410 \times 10^{-8}$	$8.250193\times 10^{-8}$
$10^{-6}$	$6.761866\times 10^{-9}$	$2.504722\times 10^{-8}$
$10^{-7}$	$2.348227\times 10^{-7}$	$2.424725\times 10^{-7}$
$10^{-8}$	$1.194055\times 10^{-6}$	$2.424645\times 10^{-6}$
$10^{-9}$	$1.763533\times 10^{-5}$	$2.424645\times 10^{-5}$
$10^{-10}$	$7.233669\times 10^{-6}$	$2.424645\times 10^{-4}$
$10^{-11}$	$1.704020\times 10^{-4}$	$2.424645\times 10^{-3}$
$10^{-12}$	$1.757870\times 10^{-2}$	$2.424645\times 10^{-2}$
$10^{-13}$	$1.493988\times 10^{-1}$	$2.424645\times 10^{-1}$
$10^{-14}$	$2.230919$	$2.424645$
$10^{-15}$	$1.875648$	$24.24645$
$10^{-16}$	$218.3926$	$242.4645$

$\Box$