Deriving Finite Difference Approximations

Univariate Case

Theory

Fix some small positive $\Delta x$ . Say we want to approximate $f^{(k)}(x)$ (i.e. the $k$ th derivative of $f$ at $x$ ), using the $d$ distinct samples $f(x + s_1 \Delta x)$ , $\ldots$ , $f(x + s_d \Delta x)$ . (Usually all $s_i$ will be integers, but this isn’t critical.) Specifically, we’d like to find coefficients $a_1, \ldots, a_d$ so that

$\sum_{n = 1}^d a_n f(x + s_n \Delta x) \approx f^{(k)}(x)$

Each sample can be expanded in a Taylor series about $x$ .

$f(x + s_i \Delta x) = \sum_{n = 0}^{\infty} \frac{(s_i \Delta x)^n}{n!} f^{(n)}(x)$

So the linear combination of samples can be written as

$\begin{aligned} \sum_{n = 1}^d a_n f(x + s_n \Delta x) &= \sum_{n = 1}^d a_n \sum_{m = 0}^{\infty} \frac{(s_n \Delta x)^m}{m!} f^{(m)}(x) \\ &= \sum_{m = 0}^{\infty} \left( \sum_{n = 1}^d a_n s_n^m \right) \frac{(\Delta x)^m}{m!} f^{(m)}(x) \end{aligned}$

The only term of the outer sum that contains $f^{(k)}(x)$ is $m = k$ . To ensure that its coefficient is one, we impose

$\sum_{n = 1}^d a_n s_n^k = \frac{k!}{(\Delta x)^k}$

In an ideal world, for all $m \neq k$ , we could impose

$\sum_{n = 1}^d a_n s_n^m = 0$

This would make our approximator perfect (i.e. not an approximation at all). But we can’t expect to satisfy an infinite number of equations with only $d$ unknowns, so we have to restrict the $m$ that we apply it to. In particular, as long as $k < d$ , we can impose this constraint for all (non-negative) $m$ such that $m \neq k$ and $m < d$ (since then we have $d$ unknown coefficients and $d$ equations). Then our approximation expands as

$\sum_{n = 1}^d a_n f(x + s_n \Delta x) = f^{(k)}(x) + R_x(\Delta x)$

where

$R_x(\Delta x) = \sum_{m = d}^{\infty} \left( \sum_{n = 1}^d a_n s_n^m \right) \frac{(\Delta x)^m}{m!} f^{(m)}(x)$

At first glance, this remainder appears to be order $d$ in $\Delta x$ . But the coefficients we solved for could have factors of $\Delta x$ that change this. Indeed, they’ll usually each have a factor of $(\Delta x)^{-k}$ . So the remainder ends up being order $d - k$ , which is at least one (recall that $k < d$ so that the system isn’t overdetermined). More often than you might expect, we get lucky and $\sum_{n = 1}^d a_n s_n^m = 0$ for some $m \geq d$ , which makes the order greater than $d - k$ . (Try using $s_1 = -1$ and $s_2 = 1$ .) To guarantee that we get a higher order approximation, we can always add more samples (thus increasing $d$ , and the number of constraints we can impose).

Implementation

The constraints that we imposed form a system of linear equations that we can write as

$\mathbf{M} \mathbf{a} = \mathbf{b}$

where $M_i^j = s_j^{i - 1}$ , $a_i$ are the unknown coefficients as defined above, and

$b_i = \begin{cases} \frac{k!}{(\Delta x)^k} &\text{if } i = k \\ 0 &\text{if } i \neq k \end{cases}$

Note that $\mathbf{M}$ is a $d \times d$ matrix, and both $\mathbf{a}$ and $\mathbf{b}$ are vectors of length $d$ .

The rows of $\mathbf{M}$ will always be linearly independent. To prove this, consider the polynomial

$g(t) = \sum_{n=1}^d \alpha_n t^{n - 1}$

As long as at least one $\alpha_n$ is nonzero, $g$ is a nonzero polynomial of order at most $d - 1$ , and thus has at most $d - 1$ distinct zeros. Thus there is no set of coefficients $\alpha_n$ such that $g(s_j) = 0$ for all $d$ samples. In other words,

$\sum_{n=1}^d \alpha_n s_j^{n - 1} \neq 0$

for at least one sample $s_j$ , as long as some $\alpha_n$ is nonzero. This is the definition of linear indedence for the rows of $\mathbf{M}$ . As such the system will always be consistent. It can be solved using an LU decomposition.

Multivariate Case

Theory

The same basic theory works in the multivariate case, the algebra just gets a little more complex. Let $f : \mathbb{R}^r \rightarrow \mathbb{R}$ be a function of $r$ variables. Given a multi-index $\alpha$ in $\mathbb{Z}_+^r$ , define

$|\alpha| = \sum_{n=1}^r |\alpha_n| = |\alpha_1| + \cdots + |\alpha_r|$ $\alpha! = \prod_{n=1}^r \alpha_n! = \alpha_1! \cdots \alpha_r!$ $\mathbf{x}^\alpha = \prod_{n=1}^r x_n^{\alpha_n} = x_1^{\alpha_1} \cdots x_r^{\alpha_r}$ $D^\alpha f = \frac{\partial^{|\alpha|}}{\partial x_1^{\alpha_1} \cdots \partial x_r^{\alpha_r}} f$

Then the multivariate Taylor expansion of a sample $f(\mathbf{x} + \mathbf{s} \odot \mathbf{\Delta x})$ about $\mathbf{x}$ (where $\mathbf{s} \odot \mathbf{\Delta x}$ denotes the element-wise or Hadamard product) is

$\begin{aligned} f(\mathbf{x} + \mathbf{s} \odot \mathbf{\Delta x}) &= \sum_{\alpha} \frac{(\mathbf{s \odot \Delta x})^\alpha}{\alpha!} D^\alpha f(\mathbf{x}) \\ &= \sum_{\alpha} \frac{\mathbf{s}^\alpha (\mathbf{\Delta x})^\alpha}{\alpha!} D^\alpha f(\mathbf{x}) \\ &= \sum_{n=0}^\infty \sum_{|\alpha| = n} \frac{\mathbf{s}^\alpha (\mathbf{\Delta x})^\alpha}{\alpha!} D^\alpha f(\mathbf{x}) \end{aligned}$

We can restrict this to a certain order by restricting the order of $\alpha$ .

So the linear combination of $d$ samples $f(\mathbf{x} + \mathbf{s_1} \odot \mathbf{\Delta x}), \ldots, f(\mathbf{x} + \mathbf{s_d} \odot \mathbf{\Delta x})$ becomes

$\begin{aligned} \sum_{n = 1}^d a_n f(\mathbf{x} + \mathbf{s_n} \odot \mathbf{\Delta x}) &= \sum_{n = 1}^d a_n \sum_{m=0}^\infty \sum_{|\alpha| = m} \frac{\mathbf{s_n}^\alpha (\mathbf{\Delta x})^\alpha}{\alpha!} D^\alpha f(\mathbf{x}) \\ &= \sum_{m = 0}^{\infty} \sum_{|\alpha| = m} \left( \sum_{n = 1}^d a_n \mathbf{s_n}^\alpha \right) \frac{(\mathbf{\Delta x})^\alpha}{\alpha!} D^\alpha f(\mathbf{x}) \end{aligned}$

And to approximate a particular mixed partial derivative $D^\kappa f(\mathbf{x})$ we impose the constraints

$\sum_{n = 1}^d a_n \mathbf{s_n}^\kappa = \frac{\kappa!}{(\mathbf{\Delta x})^\kappa}$

and, for as many $\alpha \neq \kappa$ as possible,

$\sum_{n = 1}^d a_n \mathbf{s_n}^\alpha = 0$

Finally, if instead of approximating a single mixed partial derivative, we’d like to approximate a linear combination of them, we can encode this directly into our scheme. This is useful, for example, for approximating the Laplacian of a function.

Implementation

At the end of the day, we’re still solving a system of linear equations that we can write as

$\mathbf{M} \mathbf{a} = \mathbf{b}$

where $a_i$ are the unknown coefficients. But defining $\mathbf{M}$ and $\mathbf{b}$ is a little trickier.

In particular, in the univariate case, to constrain derivatives zero through $k$ , we needed to satisfy $k + 1$ equations. But in the multivariate case, the number of mixed partial derivatives of order $k$ grows nonlinearly. In fact, the number of multi-indices $\alpha$ in $\mathbb{Z}_+^r$ of order $k$ is the binomial coefficient

$\left| \{ \alpha \text{ such that } |\alpha| = k \} \right| = \left( \frac{r + k - 1}{k} \right)$

And the number of multi-indices with order less than or equal to $k$ is

$\left| \{ \alpha \text{ such that } |\alpha| \leq k \} \right| = \left( \frac{r + k}{k} \right)$

Additionally, the multi-index exponentiation introduces a new way for the rows of $\mathbf{M}$ to be linearly dependent.

To get around these issues, we manually choose $d$ multi-indices $\beta_1, \ldots, \beta_d$ corresponding to the $d$ mixed partial derivatives $D^{\beta_1}, \ldots, D^{\beta_d}$ that we want to constrain. One of these $\beta_i$ should be $\kappa$ , certainly (the derivative we want to approximate). And the approximation is useless unless all mixed partials of order less than or equal to $\vert \kappa \vert$ are constrained to be zero, so they need to be in the mix as well.

Then

$M_i^j = \mathbf{s_j}^{\beta_i}$

and

$b_i = \begin{cases} \frac{\kappa!}{(\mathbf{\Delta x})^\kappa} &\text{if } \beta_i = \kappa \\ 0 &\text{if } \beta_i \neq \kappa \end{cases}$

If the resulting matrix $\mathbf{M}$ is singular, just swap out one constraint for another. Note that you should also pay attention to your $b_i$ : if the values are the same for two linearly dependent rows of $\mathbf{M}$ , then that constraint will still be satisfied even after it is removed (the system was underdetermined). But if the values of $b_i$ are different, then your system is inconsistent. This likely means you’re not respecting some symmetry of the problem (assuming you haven’t made a typo).

Luckily, beyond this annoyance, having $\mathbf{M}$ be singular is actually helpful for us. In the general case, we’d need $r + \vert \kappa \vert$ choose $\vert \kappa \vert$ samples in order to constrain all derivatives up to order $\vert \kappa \vert$ . But if these equations end up linearly dependent, it means they’ll all be satisfied and we only need to include a subset of them in $\mathbf{M}$ . This means we can include some higher order derivatives as well, or reduce the number of samples.

Analyzing Finite Difference Approximations

To stay general, I will only consider the multivariate case.

First, an obvious point that’s worth stating explicitly. The coefficients $a_i$ in an approximation

$\sum_{n = 1}^d a_n f(\mathbf{x} + \mathbf{s_n} \odot \mathbf{\Delta x}) \approx D^\kappa f(\mathbf{x})$

do not depend on $x$ , so we can move $x$ arbitrarily. This is essential, since it means we can apply the same approximation to every point in a discretized PDE.

Second, recall that the expansion of an approximation can be written as

$\sum_{n = 1}^d a_n f(\mathbf{x} + \mathbf{s_n} \odot \mathbf{\Delta x}) = \sum_{m = 0}^{\infty} \sum_{|\alpha| = m} \left( \sum_{n = 1}^d a_n \mathbf{s_n}^\alpha \right) \frac{(\mathbf{\Delta x})^\alpha}{\alpha!} D^\alpha f(\mathbf{x})$

So the order of the approximation can be found by finding the lowest order $\alpha \neq \kappa$ such that

$\left( \sum_{n = 1}^d a_n \mathbf{s_n}^\alpha \right) \neq 0$

However this $\alpha$ doesn’t give you the order of the error directly, since the coefficients will have powers of $\mathbf{\Delta x}$ in their denominators. So you need to divide it out to find the actual order of the error.

The leading error term is (again, consider the powers of $\mathbf{\Delta x}$ in the coefficients)

$\left( \sum_{n = 1}^d a_n \mathbf{s_n}^\alpha \right) \frac{(\mathbf{\Delta x})^\alpha}{\alpha!} D^\alpha f(\mathbf{x})$

Note that there may be multiple $\alpha$ of the order of the error; you may want to sum over them.