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. Say for example you start with the first-order derivative of some function \(f(x)\) as the 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 no limit is taken, the difference between those function values remains finite \(f(x+h) – f(x)\), therefore it is called “Finite Difference Method”.

The example above is a so-called “forward finite difference”, since one of the function values is advanced by \(h\) and the corresponding numerical finite difference scheme is named “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 allows for solving differential equations rather quickly, but might suffer from numerical instabilities (such as violation of energy conservation) and greater numerical errors on the order of \(\mathcal{O}(h)\)

### Central Difference Scheme

Therefore, in practice usually 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 = \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 slightly more elaborate, have the huge 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 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

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

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