One way of solving a differential equation is to discretize it on a lattice and approximate the derivatives using the finite difference method.

The method derives directly from the mathematical definition of a derivative. For example, if you start with the first-order derivative of some function $f(x)$ as the following limit:

$$\frac{df}{dx}(x) = f^\prime(x) = \lim_{h \to 0} \frac{f(x+h) – f(x)}{h}$$

Then $f^\prime(x)$ can be approximated by the fraction on the right-hand side, choosing an appropriately small $h$.

$$f^\prime(x) \approx \frac{1}{h} (f(x+h) – f(x)) + \mathcal{O}(h)$$

Since we don’t take the limit, the difference between those function values remains finite $f(x+h) – f(x)$, and therefore, it is called the “Finite Difference Method.”

The example above is a so-called “forward finite difference” since one of the function values is advanced by $h$. The corresponding numerical finite difference scheme is called the “Euler method.” A variation of this is the “backward finite difference.”

$$f^\prime(x) \approx \frac{1}{h} (f(x) – f(x-h)) + \mathcal{O}(h)$$

Both schemes allow for solving differential equations rather quickly but might suffer from numerical instabilities (such as violation of energy conservation) and more significant numerical errors on the order of $\mathcal{O}(h)$

## Central Difference Scheme

Therefore, in practice, slightly more advanced finite difference schemes are employed, such as the central difference scheme

$$f^\prime(x) \approx \frac{1}{\tilde{h}} (f(x+\tilde{h}/2) – f(x-\tilde{h}/2)) + \mathcal{O}(h^2)$$

or, by choosing $h = \frac{\tilde{h}}{2}$, we obtain

$$f^\prime(x) \approx \frac{1}{2h} (f(x+h) – f(x-h)) + \mathcal{O}(h^2)$$

Central differences, despite being more elaborate, have the great advantage of being symmetric in $\pm h$ and giving a more accurate approximation, with the numerical error being of order $\mathcal{O}(h^2)$.

What about second-order derivatives? We can approximate those from our first order approximation via

\begin{aligned} \frac{d^2f}{dx^2}(x) &= f^{\prime\prime}(x) = \frac{d}{dx}\frac{df}{dx} \approx \frac{d}{dx} \left( \frac{f(x+h) – f(x-h)}{2h} \right) \\ &\approx \frac{1}{2h} \left( \frac{f(x+2h) – f(x)}{2h} – \frac{f(x) – f(x-2h)}{2h} \right) \\ &= \frac{1}{(2h)^2}\left(f(x+2h) – 2f(x) + f(x-2h) \right) \end{aligned}

or, alternatively, by redefining $h$ like above

$$f^{\prime\prime}(x) = \frac{1}{h^2} \left(f(x) – 2f(x) + f(x) \right)$$

In the same way, mixed derivatives are handled:

$$\frac{\partial^2 f}{\partial x \partial y} \approx \frac{1}{2h} \left(\left(\frac{\partial f}{\partial y}\right){x+h, y} – \left(\frac{\partial f}{\partial y}\right){x-h, y}\right)$$

with

$$\left(\frac{\partial f}{\partial y}\right)_{x+h, y} \approx \frac{f(x+h, y+h) – f(x+h, y-h)}{2 h}$$

$$\left(\frac{\partial f}{\partial y}\right)_{x-h, y} \approx\frac{f(x-h, y+h) – f(x-h, y-h)}{2 h}$$

finally resulting in

\begin{aligned} \frac{\partial^2 f}{\partial x \partial y} = & \frac{1}{4h} (f(x+h, y+h) – f(x+h, y-h) \\ & – f(x-h, y+h) + f(x-h, y-h))\end{aligned}

## Affine Coordinates

Last but not least, let’s get a bit fancy and construct a central differences scheme in two dimensions for affine coordinates. These are described by a pair of covariant vectors $\vec{e}_1$, $\vec{e}_2$ forming the metric tensor

$$\hat{g}_{i,j} = \vec{e}_1 \cdot \vec{e}_2$$

which governs all distances in affine coordinates. Its inverse, the dual metric tensor $\hat{g}^{i,j} = (\hat{g}_{i,j})^{-1}$, is employed for the actual coordinate transformation

$$\begin{pmatrix} \xi^1 \ \xi^2 \end{pmatrix} = \hat{g}^{\alpha\beta} \cdot \begin{pmatrix} x \ y \end{pmatrix}$$

The first-order derivative, given by the gradient in affine coordinates, is analogous to the one in Cartesian coordinates

$$\nabla f(\xi_1, \xi_2) = \frac{1}{2 h} \begin{pmatrix} f(\xi_1+h,\xi_2) – f(\xi_1-h,\xi_2) \\ f(\xi_1,\xi_2+h) – f(\xi_1,\xi_2-h)\end{pmatrix}$$

For the Laplacian, also the dual metric tensor needs to be taken into account

\begin{aligned} \Delta f &= \nabla_i \nabla^i f = g^{i,j} \nabla_i \nabla_j f \\ &= g^{1,1} \frac{\partial^2 u}{\partial \xi_1 \partial \xi_1} +g^{1,2} \frac{\partial^2 f}{\partial \xi_1 \partial \xi_2} +g^{2,1} \frac{\partial^2 f}{\partial \xi_2 \partial \xi_1} +g^{2,2} \frac{\partial^2 f}{\partial \xi_2 \partial \xi_2} \end{aligned}

resulting in

\begin{aligned} h^2 \cdot \Delta f &=g^{1,1} (f(\xi_1+h, \xi_2) – 2 f(\xi_1, \xi_2) + f(\xi_1-h, \xi_2)) \\ & + g^{2,2} (f(\xi_1, \xi_2+h) – 2 f(\xi_1, \xi_2) + f(\xi_1, \xi_2-h)) \\ & + \frac{g^{1,2}}{2} (f(\xi_1+h, \xi_2+h) – f(\xi_1+h, \xi_2-h) \\ & – f(\xi_1-h, \xi_2+h) + f(\xi_1-h, \xi_2-h))\end{aligned}

where we have used $g^{1,2} = g^{2,1}$

## Ressources

www.iue.tuwien.ac.at/phd/heinzl/node27.html

www.holoborodko.com/pavel/numerical-methods/numerical-derivative/central-differences/#comment-18109

www.mathematik.uni-dortmund.de/~kuzmin/cfdintro/lecture4.pdf

Published in Numerics