Numerical Differentiation:

The basic idea behind numerical differentiation is very simple. If given the set of values $(x_i, f_i)$ i=0,1,...,n, we determine the interpolating polynomial $p(x)$ through these points. We then differentiate this polynomial to obtain $p'(x)$ whose values for any $x$ is taken as an approximation to $f'(s)$. Let us very briefly describe this interpolating polynomial.
Let $x_i,\,\, i=0,1,...,n$ be $n+1$ distinct points on an interval I and let $f(x)$ be a real valued function which takes on the values $f(x_i)=f_i,i=0,1,..,n$, at these n+1 points. To construct a polynomial of degree not exceeding n which passes through the n+1 points $(x_i, f_i),\,\,i=0,1,..n,$ we use the method of Lagrange. We begin by expressing the desired polynomial as

$$p(x)=L_0(x)f_0+L_1(x)f_1+...+L_nf_n=\sum\limits_{k=0}^nL_K(x)f_k$$

where each $L_K(x),\,k=0,1,...n$ is a polynomial of degree not exceeding n. The interpolation condition $p(x_i)=f_i$ will be satisfied if the $L_K(x)$ satisfy:

$$L_K(x_i)=1\qquad K=i$$

$$\qquad\quad=0 \qquad K \neq i$$

It is easy to verify that the function defined by

$$L_K(x)=\frac{(x-x_0)(x-x_1)...(x-x_{K-1})(x-x_{K+1})...(x-x_n)}{(x_k-x_0)(x_K-x_1)...(x_K-x_{K-1})(x_K-x_{K+1})...(x_K-x_n)}$$

have this property.
The error in approximating $f(x)$ by such a polynomial $p(x)$ is given by

$$f(x)-p(x)=\frac{\psi (x)f^{(n+1)}_{\xi}}{(n+1)!}$$

where and $\psi (x)=(x-x_0)(x-x_1)...(x-x_n)$ and $\psi$ is some point on the interval containing $x_i,i=0,1,...,n$ Then the error in the derivative $p'(x)$ at a tabular point $x=x_i$ is given by

$$f'(x_i)-p'(x_i)=\frac{\psi'(x_i)f^{(n+1)}_{(\xi_i)}}{(n+1)!}$$

Where $\xi_i$ is a point in the interval containing the points $x_i\, ,i=0(1)n$ and

$$\psi'(x_i)=\prod\limits_{j=0\,j\neq i}^{n}(x_i-x_j)$$

If we assume that the interpolating points $x_{i-1} i=0(1)n$ are equally spaced with spacing h, and put $S=\frac{x-x_0}{h}$ we can approximate $f(x)$ by Newton's forward difference formula given by

$$f(x)\simeq P(s)=f_0+S\Delta f_0+\left(% \begin{array}{c} S... ...gin{array}{c} S \\ n \\ \end{array}% \right)\Delta^nf_0$$

Where $\Delta$ is the forward difference operation definition

$$\Delta f(x_0)=f(x_0+h)-f(x_0)$$

$$\Delta^kf(x_0)=\Delta^{k-1}f(x_0+h)-\Delta^{k-1}f(x_0)\qquad k=2,3,...$$

and $$\left(% \begin{array}{c} S \\ n \\ \end{array}% \right)$$ is the binomial coefficient given by

$$\left(% \begin{array}{c} S \\ n \\ \end{array}% \right)=\frac{S(S-1)(S-2)...(S-n+1)}{n!}$$

The error is $$f(x)-p(s)=\left(% \begin{array}{c} S \\ n+1 \\ \end{array}% \right)h^{(n+1)}f^{(n+1)}(\xi)$$.
Differentiating (1) w. r. h. x and noting that $\frac{dP}{dx}=\frac{1}{h}\frac{dP}{dS}$
We obtain

$$f'(x)\simeq\frac{1}{h}\frac{dP}{dS}=\frac{1}{h}[\Delta f_0+\f... ...in{array}{c} S \\ n \\ \end{array}% \right)\Delta^nf_0]$$

If we now set $x=x_0, \,($ie$\,\, S =0)$ we obtain Approximation formulas for $f'(x_0)$ for different values of n. For example, for $n=1$, we have

$$f'(x_0)\simeq\frac{\Delta f_0}{h}=\frac{f_1-f_0}{h}$$

with the error of this formula given by

$$E=\frac{\psi' (x_0)f''(\xi)}{2}=-\frac{h}{2}f''(\xi)$$

Taking now n=2,

$$f'(x_0)\approx\frac{1}{h}(\Delta f_0-\frac{1}{2}\Delta^2f_0)$$

$$=\frac{1}{2h}(-3f_0+4f_1-f_2)$$

with the error given by

$$E=\frac{\psi'(x_0)}{3!}f''(\xi)=\frac{(x_0-x_1)(x_0-x_2)}{6}f''(\xi)=\frac{h^2}{3}f'''(\xi)$$

Formulas for approximation higher derivation of $f(x)$ can be obtained in a similar manner. Differentiating (1) twice and setting $\xi=0$, we obtain

$$f''(x_0)\cong\frac{\Delta^2f_0}{h^2}=\frac{f_2-2f_1++f_0}{h^2}$$

with the error given by $E=-hf'''(\xi)$
Another formula for the derivatives, using central difference, is given by

$$f'(x_0)\simeq \frac{f_1-f_{-1}}{2h}$$

where $f_1=f(x_0+h); f_{-1}=f(x_0-h)$
Using Taylor's expansion about the $x_0$, we get

$$f_1=f_0+hf'_0+\frac{h^2}{2}f''_0+\frac{h^3}{6}f'''(\xi_1)\qquad x_0 < \xi_1 < x_1$$

$$f_{-1}=f_0-hf'_0+\frac{h^2}{2}f''_0-\frac{h^3}{6}f'''(\xi_2)\qquad x_{-1} < \xi_2 < x_0$$

This gives

$$\frac{f_1-f_{-1}}{2h}=f'(x_0)+\frac{h^2}{12}[f'''(\xi_1)+f'''(\xi_2)]$$

If $f''(x)$ is continuous on $[x_{-1},x_1]$, we have

$$f'(x_0)=\frac{f_1-f_{-1}}{2h}-\frac{h^2}{6}f'''(\xi)\,, \qquad x_{-1}<\xi < x_1$$

as thus the error in the formula is

$$E=-\frac{h^2}{6}f'''(\xi)$$

The simplest formula $fx f''(x_0)$ based on central difference is

$$f''(x_0)\simeq \frac{f_1-2f_0+f_{-1}}{h^2}$$    with error $$E=-\frac{h^2}{12}f''(\xi)$$