# Finite Difference Method

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

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

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