Numerical Differentiation

Introduction

Often we need to approximate the derivative of a function when we cannot obtain it analytically. Here we discuss several ways to do this numerically.

Taylor Series

Let's consider the situation where we have samples of a function $f(x)$ at discrete points $\begin{array}{cccc}x_1 & x_2 & \cdots & x_n\end{array}$ seperated by spacing $h$ as depicted in the following figure:

Consider a Taylor series expansions about some arbitrary point $x_i$. Since $x_{i+1}-x_i = x_i-x_{i-1} = h$ we can write these as follows:

Approximation location	Taylor Series Expansion about $x_i$
$x_{i+1}$	$f(x_{i+1}) = f(x_i) + f^\prime(x_i) h + \tfrac{1}{2}f^{\prime\prime}(x_i) h^2 + \tfrac{1}{6}f^{\prime\prime\prime}(x_i) h^3 + \tfrac{1}{24}f^{(4)}(x_i) h^4 + \cdots$
$x_{i-1}$	$f(x_{i-1}) = f(x_i) - f^\prime(x_i)h + \tfrac{1}{2}f^{\prime\prime}(x_i)h^2 - \tfrac{1}{6}f^{\prime\prime\prime}(x_i)h^3 + \tfrac{1}{24}f^{(4)}(x_i)h^4 - \cdots$

If we subtract the Taylor series expansion at $x_{i-1}$ from the one at $x_{i-1}$, we find

$f(x_{i+1})-f(x_{i-1}) = 2f^{\prime}(x_i)h + \tfrac{1}{3} f^{\prime\prime\prime}(x_i) h^3 + \cdots$

Now we solve this for $f^\prime(x_i)$ to find

$f^{\prime}(x_i) = \frac{f(x_{i+1})-f(x_{i-1})}{2h} - \tfrac{1}{6} f^{\prime\prime\prime}(x_i)h^2 - \cdots$

Now if $h$ is small, then the second term (with the $h^2$ in it) is small and we can approximate the derivative as

$f^\prime(x_i) = \frac{f(x_{i+1})-f(x_{i-1})}{2h}$

We call this a second order approximation to $f^\prime(x)$ because when we truncated the series approximation to $f^\prime(x_i)$ the largest term there was of the order of $h^2$.

Note that we now have a way to approximate the derivative of a function if we have the function's values at two locations.

On a uniform mesh, we can use this technique to generate a variety of approximations to derivatives, as summarized in the following table:

Formula for $\left. \frac{\mathrm{d} f}{\mathrm{d}x} \right|_{i}$ Order $\frac{f(x_{i+1})-f(x_{i})}{h}$ $\mathcal{O}\left( h \right)$ $\frac{f(x_{i})-f(x_{i-1})}{h}$ $\mathcal{O}\left( h \right)$ $\frac{f(x_{i+1})-f(x_{i-1})}{2 h}$ $\mathcal{O}\left(h^2 \right)$ $\frac{f(x_{i-2}) - 8 f(x_{i-1}) + 8 f(x_{i+1}) - f(x_{i+1})}{12 h}$ $\mathcal{O}\left( h^4 \right)$

Error in approximation of $\frac{df}{dx}$ as a function of the size of the interval $h$. Click here to download the matlab file that produced this plot.

As can be seen in the figure above, higher order approximations result in significantly lower error for a given spacing $h$. Note that for $h<10^{-3}$ the fourth-order approximation is contaminated by roundoff error. The same would happen for the other derivative approximations, but at smaller $h$.

Lagrange Polynomials

Lagrange polynomials, which are commonly used for interpolation, can also be used for differentiation. The formula is

$f^{\prime}(x) = \sum_{k=0}^n y_k L_{k}^{\prime}(x),$

where $L_{k}^{\prime}(x)$ is given as

$L_{k}^{\prime}(x) = \left[ \sum_{{j=0}\atop{j\ne k}}^{n} \left( \prod_{ {i=0}\atop{i \ne j,k} }^n (x-x_i) \right) \right] \left[ \prod_{{i=0}\atop{i\ne k}}^{n} (x_k-x_i) \right]^{-1}.$

Here $n$ is the order of the polynomial and we require $n_p=n+1$ points to form the Lagrange polynomial.