## Jun 20, 2014

### Numerical Differentiation

Let $f(x)$ be a smooth function. We want to take its derivative at some point $x$, numerically. In this post, we assume the function is sampled on an equal length mesh. $h$ is the step, $0 \le h \ll 1$. We consider the leading order and higher order methods. Here by order we mean the order of residual error, $\mathcal O(h^n)$.

There are three standard way: the forward difference, the backward difference and the central (or symmetric) difference. They are defined as following: $f'(x) \simeq \frac{f(x+h)-f(x)}{h} + \mathcal O (h), \\ f'(x) \simeq \frac{f(x)-f(x-h)}{h} + \mathcal O (h), \\ f'(x) \simeq \frac{f(x+h)-f(x-h)}{2h} + \mathcal O (h^2)$ The numerical accuracy can be seen from the Taylor expansion: $f(x+h) = f(x) + h f'(x) + \frac{h^2}{2}f''(x) + \frac{h^3}{3!}f'''(x) + ... \\ f(x-h) = f(x) - h f'(x) + \frac{h^2}{2}f''(x) - \frac{h^3}{3!}f'''(x) + ....$ This method, used recursively, also produces higher order numerical differentiation. Central difference for example: $f^{(n)}(x) \simeq \frac{f^{(n-1)}(x+h) - f^{(n-1)}(x-h)}{2h}$ Therefore [theorem 1], $f^{(n)}(x) = \frac{1}{h^n}\sum_{k = 0}^n {n \choose k} (-1)^k f(x+(n/2-k)h) + \frac{n h^2}{4!}f^{(n+2)}(\xi);$ The accuracy can be proven by substituting in the Taylor expansion. Just note (Stirling number cf. The Stirling number S(n,k) is the number of ways to partition a set of n objects into k groups.) $S(m,n) \equiv \frac{1}{n!}\sum_{k=0}^n {n \choose k} (-1)^{n-k} k^m = \left\{ \begin{array}{ll} 0, & m<n \\ 1, &m=n \\ n(n+1)/2, &m=n+1\\ S(n-1,m-1) + m S(n-1,m), & \\ \end{array}\right.$ When $n = 1 \bmod 2$, the meshes are on half integers. One can use $h \to 2h$ to recover the integral grid. The series sometimes is defined as an operator known as the central difference operator (Abramowitz and Stegun, $\S$25.1, "Differences"): $\delta_n^{2k} f = \sum_{j=0}^{2k} (-1)^j { 2k \choose j} f_{n+k-j}, \\ \delta_{n+\frac{1}{2}}^{2k+1} f = \sum_{j=0}^{2k+1} (-1)^j { 2k+1 \choose j} f_{n+k+1-j}.$ Then our theorem 1 can be written as: $\frac{d^n}{d x^n} f(x) = \frac{\delta_0^{n}}{\Delta x^n} f(x) + \mathcal O(\Delta x^2)$.

#### High Order Finite Differences

In the previous section, we have discussed approximate the differentiation numerically up to leading order $\mathcal O(h)$. The higher order numerical differentiations are rarely used but nevertheless helpful to know. Inspired by the central difference series, in general, we can approximate the $n$-th order derivative by series $f^{(n)}(x) = \sum_k c_k^{(n,r)} f(x + k h) + \mathcal O(h^r)$ Obviously, the larger $n$ and/or $r$ is, the more terms we need. We can use the Lagrange polynomial to approximate the function $f(x)$ up to $f^{(n+r)}(\xi)$ and then take the $n$-th order derivative of the Lagrange polynomial. Then $f^{(n)} (x) = \sum_k l_k^{(n)}(x) f(x+kh),$ where $l_i(x) = \prod_{j\ne i} \frac{x-x_j}{x_i-x_j}$.

example: the finite-difference approximation of the first order derivative ($f_n = f(x+nh)/h$): $f'(x) \simeq \sum_{k=1}^{n} \frac{(-1)^{k-1} (n)_k}{(n+1)^{(k)}}\frac{f(x+kh) - f(x-kh)}{kh},$ where $(x)_k = x(x-1)\cdots(x-k+1)$ is the falling factorial and $x^{(k)} = x(x+1)\cdots(x+k-1)$ is the rising factorial.
second order derivatives ($f_n = f(x+nh)/h^2$): $f''(x) \simeq 2\sum_{k=1}^n \frac{(-1)^{k-1} (n)_k}{(n+1)^{(k)}} \frac{ f(x+kh)+f(x-kh) }{k^2h^2} - \frac{2H_{n,2}}{ h^2} f(x),$ where $H_{n,2} = \sum_{k=1}^n \frac{1}{k^2}$ is called harmonic number.
third order derivatives ($f_n = f(x+nh)/h^3$):

In Mathematica, function InterpolationPolynomial[] provides the Lagrange polynomial. The coefficients can be computed as following:
L[K_] := InterpolatingPolynomial[
Table[{k, Subscript[f, k]}, {k, -K, K}], x]

fd[K_, n_] := Collect[#, Variables[#]] & @(D[L[K], {x, n}] /. x -> 0)