Site Overlay

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

Ben's Planet is built upon CC-BY-SA
Scroll Up