A previous post used the typical series-expansion-plus-linear-algebra approach for finding finite-difference formulas to derive approximations to the first derivative of any desired order of accuracy. If you’ve ever used that method yourself, you probably know how tedious it is, and that post doesn’t make it look much less so. On top of that, for all that work we only got approximations for central differences. If you wanted approximations with more points on one side than another, you would have to go through the whole process again.
The bright side of the method in that post was that no linear system was solved, which is the annoying part of finding these formulas by hand. We avoided that step by noticing that the solution to the system could be read off from the Lagrange basis polynomials. So one might suspect that focusing more on the Lagrange basis from the outset might make deriving finite-difference formulas even easier, and that is the subject of this post.
Foundational Ideas
For a given set of points $\{x_i\}_{i=-a}^{i=b}$, the Lagrange basis polynomial $L_j(x)$ associated with $x_j$ is the polynomial which equals one at $x_j$ and zero at all the other $x_i$.
$$\displaystyle L_j(x_i) = \delta_{ij} = \begin{cases} 1 & i = j \\ 0 & i \neq j \end{cases} $$
This definition makes it very easy to write the polynomial $P(x)$ that interpolates values $f_i$ defined at the $x_i$:
$$\displaystyle P(x) = \sum_j f_jL_j(x) $$
$L_j(x)$ can be written explicitly in terms of the node locations $x_j$:
$$\displaystyle L_j(x) = \frac{\prod_{k\neq j} (x – x_j)}{\prod_{k\neq j} (x_k – x_j)} = \frac{\prod_{k\neq j} (x – x_j)}{D_j}$$
The finite-difference formulas we seek are those which are exact for all polynomials of sufficiently low degree, in particular they will be exact for the interpolant $P(x)$. So for some arbitrary unknown function $f(x)$, we will have:
$$\displaystyle \begin{aligned} f(x) &\approx P(x) = \sum_j f_j L_j(x) \\ \frac{d^n f}{dx^n}\Bigg|_{x=0} &\approx \sum_j f_j \left( \frac{d^n L_j}{dx^n}\Bigg|_{x=0} \right) = \sum_j c_jf_j \end{aligned}$$
Thus the coefficients $c_j$ we want are the $n$-th derivatives of the Lagrange polynomials evaluated at $x=0$, which we can evaluate directly. For concreteness and simplicity, we will consider points $x_i = ih$ where $i = -a, -a+1, \cdots, 0, \cdots, b-1, b$ where $a>0$ and $b>0$ and $h$ is a fixed gap size. That is, we will assume that $x=0$ is within the interval spanned by the $x_i$.
Differentiating the $L_j$
Since the denominator $D_j$ can be factored out, let’s first get an expression for it:
$$\displaystyle D_j = \prod_{\substack{k=-a, \\ k\neq j}}^{k=b} (kh – jh) = h^{a+b}\prod_{\substack{k=-a, \\ k\neq j}}^{k=b} (k – j)$$
Since the factors are $(k-j)$ and $j$ is fixed for a given $L_j$, we can re-index the product replacing $k$ by $v = k-j$ and get:
$$\displaystyle \begin{aligned} D_j &= h^{a+b}\prod_{\substack{v=-a-j, \\ v\neq 0}}^{v=b-j} v = h^{a+b}\left( \prod_{v=-a-j}^{v=-1} v \right)\left( \prod_{v=1}^{v=b-j} v \right) \\ &= h^{a+b}(-1)^{a+j}(a+j)!(b-j)! \end{aligned}$$
So now we need the $n$-th derivative of the numerators:
$$\displaystyle A^n_j (x)= \frac{d^n}{dx^n} \prod_{k\neq j} (x – x_j) $$
Following the product rule, each of the $n$ derivatives splits products into sums of products each missing one factor from the first product:
$$\displaystyle \frac{d}{dx}\prod_{k\in S} (x-x_k) = \sum_{v\in S}\prod_{\substack{k\in S \\ k\neq v}} (x-x_k) $$
Applying this process $n$ times:
$$\displaystyle A_j^n (x) = \sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{j\} \\ |S|=n}} \left( n! \prod_{\substack{k\notin S \\ k \neq j}} (x-x_k) \right) $$
The $n!$ factor appears because there are $n!$ ways in which the $n$ derivatives could be applied to produce the same factor in the final sum. For example, the the first of the $n$ derivatives could be applied to the $x-x_1$ factor, the second derivative to the $x-x_{-3}$ factor, and so on. Or, the first derivative could be applied to $x-x_{-3}$ and the second to $x-x_1$, and so on. Each of these orderings produces a copy of the same term, and there are $n!$ orderings.
Now we obtain the coefficients by evaluating $A_j^n(x)$ at $x=0$:
$$\displaystyle \begin{aligned} A_j^n(0) &= n!\sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{j\} \\ |S|=n}} \prod_{\substack{k\notin S \\ k \neq j}} (-kh) \\ &= n!(-h)^{a+b-n} \sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{j\} \\ |S|=n}} \prod_{\substack{k\notin S \\ k \neq j}} (k)\end{aligned}$$
Examining the products being summed in this last expression, If $j\neq 0$ and $0 \notin S$, then one of the factors in the product will be zero making the whole product zero, so that it does not contribute to the sum. We have two cases depending on whether $j=0$:
$$\displaystyle \sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{j\} \\ |S|=n}} \prod_{\substack{k\notin S \\ k \neq j}} (k) = \begin{cases} j=0 & \sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{j\} \\ |S|=n}} \frac{a!b!(-1)^{a}}{\prod_{k\in S} k} \\ j \neq 0 & \sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{j,0\} \\ |S|=n-1}} \frac{a!b!(-1)^a}{j\prod_{k\in S} k} \end{cases} $$
Here we have made use of the fact that:
$$\displaystyle \prod_{k=-a}^{k=b} k = (-1)^a a!b! $$ and divide that product by the values $k$ present in each subset $S$. Putting all of this together, we obtain for $c_0$
$$\displaystyle \begin{aligned} c_0 &= \frac{A_0^n(0)}{D_0} = \frac{n!h^{a+b-n}}{h^{a+b}(-1)^a a! b!} \times a!b!(-1)^a \sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{0\} \\ |S|=n}} \frac{1}{\prod_{k\in S} k} \\ &= \frac{n!}{h^n} \sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{0\} \\ |S|=n}} \frac{1}{\prod_{k\in S} k} \end{aligned}$$
And for $j\neq 0$:
$$\displaystyle c_j = \frac{A_j^n(0)}{D_j} = \frac{(-1)^{a+b+n+j}a!b!n!}{(a+j)!(b-j)!jh^n} \sum_{\substack{S\subseteq \{-a,\cdots,b\}\setminus\{j,0\} \\ |S|=n-1}} \frac{1}{\prod_{k\in S} k} $$
An Example
Let’s compute the 5-point central difference for the second derivative. So $n = 2, a=2, b = 2$
$$\displaystyle c_{-2} = \frac{2!\cdot 2!\cdot 2!}{-2h^2\cdot 0!\cdot 4!} \left( \frac{1}{-1} + \frac{1}{1} + \frac{1}{2} \right) = \frac{-1}{12h^2}$$
$$\displaystyle c_{-1} = \frac{-2!\cdot 2!\cdot 2!}{-h^21!\cdot 3!} \left( \frac{1}{-2} + \frac{1}{1} + \frac{1}{2} \right) = \frac{4}{3h^2}$$
$$\displaystyle \begin{aligned} c_0 &= \frac{2}{h^2} \left( \frac{1}{(-2)(-1)} + \frac{1}{(-2)(1)} + \frac{1}{(-2)(2)} + \frac{1}{(-1)(1)} + \frac{1}{(-1)(2)} + \frac{1}{(1)(2)} \right) \ &= \frac{-5}{2h^2} \end{aligned}$$
$$\displaystyle c_{1} = \frac{-2!\cdot 2!\cdot 2!}{h^2 3!\cdot 1!} \left( \frac{1}{-2} + \frac{1}{-1} + \frac{1}{2} \right) = \frac{4}{3h^2}$$
$$\displaystyle c_{2} = \frac{2!\cdot 2!\cdot 2!}{2h^2\cdot 4!\cdot 0!} \left( \frac{1}{-2} + \frac{1}{-1} + \frac{1}{1} \right) = \frac{-1}{12h^2}$$
Which agrees with what one finds in the table here.
Wrapping Up
We have obtained a formula that produces the coefficients in a finite difference formula over any set of equally spaced points for derivatives of arbitrary order. The key was to use the Lagrange form of the interpolating polynomial to express the coefficients as derivatives of the Lagrange basis polynomials at $x=0$, which we then were able to compute. Hooray!