Introduction: The Finite Differences Method (FDM)

Before progressing towards finite element modeling, we will learn about the Finite Difference Method (FDM), which is a somewhat easier method to solve partial differential equations. We will do so by looking at how heat conduction “works”.

Background on conductive heat transport

Temperature is one of the key parameters that control and affect the geological processes that we are exploring during this class. Mantle melting, rock rheology, and metamorphism are all functions of temperature. It is therefore essential that we understand how rocks exchange heat and change their temperature.

One fundamental equation in the analysis of heat transfer is Fourier’s law for heat flow:

(1)\[\begin{split}\vec{q} = -k \nabla T = -k \begin{bmatrix} \frac{\partial T}{\partial x} \\ \frac{\partial T}{\partial y} \\ \frac{\partial T}{\partial z} \end{bmatrix}\end{split}\]

It states that heat flow is directed from high to low temperatures (that’s where the minus sign comes from) and is proportional to the geothermal gradient. The proportionality constant, k, is the thermal conductivity which has units of W/m/K and describes how well a rock transfers heat. k is a typically a complex function of rock type, porosity, and temperature yet is often simplified to a constant value.

In most cases, we are interested in temperature and not heat flow so that we would like to have an equation that describes the temperature evolution in a rock. We can get this by starting from the convervations equation of internal energy.

Fig. 1 Derivation of the energy equation. Change in internal energy is related to changes in heat fluxes into and out of the box (and a lot of other terms (e.g. advection) that are neglected here).

For simple incompressible cases, changes in internal energy can be expressed as changes in temperature times density and specific heat. The change in internal energy with time can now be written as \(\rho c_p \frac{\partial}{\partial t} T \Delta x \Delta y \Delta z\) (\(\rho\) is density, \(c_p\) specific heat, and \(T\) temperature), has units of (\(J/s\)), and must be equal to the difference between the heat flow into the box \(q_{in} \Delta y \Delta z \left( \frac{W}{m K}\frac{K}{m}m^2 = \frac{J}{s}\right)\) and the heat flow out of the box \(\left( q_{in} + \frac{\partial q_{in}}{\partial x} \Delta x\right) \Delta y \Delta z\) (the y and z directions are done in the same way). With these considerations, we can write a conservation equation for energy in 1D:

(2)\[\begin{split}\begin{align} \begin{split} \rho c_p \frac{\partial T}{\partial t} = -\frac{\partial q_{in}}{\partial x} = \frac{\partial}{\partial x}k\frac{\partial T}{\partial x}\\ \frac{\partial T}{\partial t} = \frac{k}{\rho c_p} \frac{\partial^2 T}{\partial x^2} = \kappa \frac{\partial^2 T}{\partial x^2}\\ \frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2} \end{split} \end{align}\end{split}\]

equation (2) is called the heat transfer or heat diffusion equation and is one of the most fundamental equations in Earth Sciences. If the thermal conductivity is constant, we can define a thermal diffusivity, \(\kappa\), and write the simpler second form of the equation.

Note that in 3D, the changes in heat flow are expressed as a divergence:

(3)\[\rho c_p \frac{\partial T}{\partial t} = -\nabla \cdot (\vec{q}) = \frac{\partial}{\partial x}k\frac{\partial T}{\partial x} + \frac{\partial}{\partial y}k\frac{\partial T}{\partial y} + \frac{\partial}{\partial z}k\frac{\partial T}{\partial z}\]

Or in vector notation:

(4)\[\begin{split}\rho c_p \frac{\partial T}{\partial t} = -\nabla \cdot \left(\vec{q}\right) = \begin{bmatrix} \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \end{bmatrix} k \begin{bmatrix} \frac{\partial T}{\partial x}\\ \frac{\partial T}{\partial y}\\ \frac{\partial T}{\partial z} \end{bmatrix}\end{split}\]

Finite Differences discretization

equation (2) is a partial differential equation that describes the evolution of temperature. There are two fundamentally different ways to solve it: 1) analytically or 2) numerically. Analytical solutions have the virtue that they are exact but it is often not possible to find one for complex systems. Numerical solutions are always approximations but can be found also for very complex systems. We will first use one numerical technique called finite differences. To learn how partial differential equations are solved using finite differences, we have to go back to the definition of a derivative:

(5)\[\frac{\partial T}{\partial x} = \lim_{\Delta x \rightarrow 0} \frac{T(x + \Delta x) - T(x)}{\Delta x}\]

In our case, the above derivative describes the change in temperature with x (our space coordinate). In numerics, we always deal with discrete functions (as opposed to continuous functions), which only exist at grid points Fig. 2. We can therefore translate the above equation into computer readable form:

(6)\[\frac{\partial T}{\partial x} \approx \frac{T(i + 1) - T(i)}{x(i+1) - x(i)} = \frac{T(i + 1) - T(i)}{\Delta x}\]

\(T(i)\) is the temperature at a grid point \(i\) and \(\Delta x\) is the grid spacing between two grid points. Using this definition of a derivative, we can now compute the heat flow from the calculated temperature solution:

(7)\[q_x = -k \frac{\partial T}{\partial x} \approx -k \left(\frac{T(i+1)-T(i)}{x(i+1)-x(i)}\right)\]

This form is called the finite differences form and is a first step towards solving partial differential equations numerically.

Note that it matters in which direction you count. Usually it makes life much easier if indices and physical coordinates point in the same direction, e.g. x coordinate and index increase to the right.

We have learned how we can compute derivatives numerically. The next step is to solve the heat conduction equation (equation (2)) completely numerically. We are interested in the temperature evolution as a function of time and space \(T(x, t)\), which satisfies equation (2), given an initial temperature distribution. We know already from the heat flow example how to calculate first derivatives (forward differencing):

(8)\[\frac{\partial T}{\partial t} = \frac{T_i^{n+1} - T_i^n}{\Delta t}\]

The index \(n\) corresponds to the time step and the index \(i\) to the grid point (x-coordinate, Fig. 2). Next, we need to know how to write second derivatives. A second derivative is just a derivative of a derivate. So we can write (central differencing):

(9)\[\kappa\frac{\partial^2 T}{\partial x^2} \approx \kappa \frac{\frac{T_{i+1}^n - T_i^n}{\Delta x}-\frac{T_{i}^n - T_{i-1}^n}{\Delta x}}{\Delta x} = \kappa \frac{T_{i+1}^n-2T_{i}^n+T_{i-1}^n}{\Delta x^2}\]

If we combine equation equation (8) and equation (9) we get:

(10)\[\frac{T_i^{n+1}-T_i^n}{\Delta t} = \kappa \left(\frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{{\Delta x^2}}\right)\]

Tip

Notice how we have conveniently used the time index \(n\) for the temperatures in the spatial derivatives. This results in the explicit form of the final discretized equation. The implicit form, which we will learn about later, would used the unknown new (time index \(n+1\)) temperatures for the spatial derivatives, which requires solving a system of equations.

The last step is a rearrangement of the discretized equation, so that all known quantities (i.e. temperature at time \(n\)) are on the right-hand side and the unknown quantities on the left-hand side (properties at \(n+1\)). This results in:

(11)\[T_i^{n+1} = \frac{\kappa \Delta t}{\Delta x^2} \left( T_{i+1}^n - 2 T_i^n + T_{i-1}^n\right) + T_i^n\]

We have now translated the heat conduction equation equation (2) into a computer readable finite differences form.

Appendix

Taylor-series expansions

Finite difference approximations can be derived through the use of Taylor-series expansions. Suppose we have a function \(f(x)\), which is continuous and differentiable over the range of interest. Let’s also assume we know the value \(f(x_0)\) and all the derivatives at \(x = x_0\). The forward Taylor-series expansion for \(f(x_0+\Delta x)\) about \(x_0\) gives

(12)\[f(x_0+\Delta x) = f(x_0)+ \frac{\partial f(x_0)}{\partial x} \Delta x + \frac{\partial^2 f(x_0)}{\partial x^2} \frac{(\Delta x)^2}{2!} + \frac{\partial^3 f(x_0)}{\partial x^3} \frac{(\Delta x)^3}{3!} + \frac{\partial^n f(x_0)}{\partial x^n} \frac{(\Delta x)^n}{n!} + O(\Delta x)^{n+1}\]

We can compute the first derivative by rearranging equation ref{eqs:Taylor_series_expansion}

(13)\[\frac{\partial f(x_0)}{\partial x} = \frac{f(x_0+\Delta x)− f(x_0)}{\Delta x} − \frac{\partial^2 f(x_0)}{\partial x^2} \frac{(\Delta x)}{2!} − \frac{\partial^3 f(x_0)}{\partial x^3} \frac{(\Delta x)^2}{3!} ...\]

This can also be written in discretized notation as:

(14)\[\frac{\partial f(x_i)}{\partial x} = \frac{f_{i+1}−f_i}{\Delta x} + O(\Delta x)\]

here \(O(\Delta x)\) is called the truncation error, which means that if the distance \(\Delta x\) is made smaller and smaller, the (numerical approximation) error decreases with \(\Delta x\). This derivative is also called first order accurate.

We can also expand the Taylor series backward

(15)\[f(x_0−\Delta x) = f(x_0)− \frac{\partial f(x_0)}{\partial x} \Delta x+ \frac{\partial^2 f(x_0)}{\partial x^2} \frac{(\Delta x)^2}{2!} − \frac{\partial^3 f(x_0)}{\partial x^3} \frac{(\Delta x)^3}{3!} + \cdots\]

In this case, the first (backward) derivative can be written as

(16)\[\begin{split}\frac{\partial f (x0)}{\partial x} = \frac{f(x_0)− f(x_0−\Delta x)}{\Delta x} + \frac{\partial^2 f(x_0)}{\partial x^2} \frac{(\Delta x)}{2!} − \frac{\partial^3 f(x_0)}{\partial x^3} \frac{(\Delta x)^2}{3!} \cdots\\ \frac{\partial f(x_i)}{\partial x} = \frac{f_i− f_{i−1}}{\Delta x} + O(\Delta x)\end{split}\]

By adding equation (13) and equation (16) and dividing by two, a second order accurate first order derivative is obtained

(17)\[\frac{\partial f(x_i)}{\partial x} = \frac{f_{i+1} − f_{i−1}}{2\Delta x} + O(\Delta x)^2\]

By adding equations equation (13) and equation (15) an approximation of the second derivative is obtained

(18)\[\frac{\partial f^2(x_i)}{\partial x^2} = \frac{f_{i+1}−2 f_i+ f_{i−1}}{(\Delta x)^2} +O(\Delta x)^2\]

With this approach we can basically derive all possible finite difference approximations. A different way to derive the second derivative is by computing the first derivative at \(i+\frac{1}{2}\) and at \(i-\frac{1}{2}\) and computing the second derivative at \(i\) by using those two first derivatives:

(19)\[\begin{split}\begin{align} \begin{split} \frac{\partial f(x_{i+1/2})}{\partial x} = \frac{f_{i+1}−f_i}{x_{i+1}−x_i}\\ \frac{\partial f(x_{i−1/2})}{\partial x} = \frac{f_i−f_{i−1}}{x_i−x_{i−1}}\\ \frac{\partial f^2(x_i)}{\partial x^2} = \frac{\frac{\partial f(x_{i+1/2})}{\partial x} − \frac{\partial f(x_{i−1/2})}{\partial x}}{x_{i+1/2}−x_{i−1/2}} = \frac{ \frac{f_{i+1}−f_i}{x_{i+1}−x_i} − \frac{f_i−f_{i−1}}{x_i−x_{i−1}}}{0.5(x_{i+1}−x_{i−1})} \end{split} \end{align}\end{split}\]

Similarly we can derive higher order derivatives. Note that the highest order derivative that usually occurs in geodynamics is the \(4^th\)-order derivative.

Finite difference approximations

The following equations are common finite difference approximations of derivatives. If you, in the future, need to write a finite difference approximation, come back here.

Left-sided first derivative, first order

(20)\[\left| \frac{\partial u}{\partial x}\right|_{i−1/2} = \frac{u_i−u_{i−1}}{\Delta x} +O(\Delta x)\]

Right-sided first derivative, first order

(21)\[\left| \frac{\partial u}{\partial x}\right|_{i+1/2} = \frac{u_{i+1}−u_i}{\Delta x} +O(\Delta x)\]

Central first derivative, second order

(22)\[\left| \frac{\partial u}{\partial x}\right|_i = \frac{u_{i+1}−u_{i−1}}{2\Delta x} +O(\Delta x)^2\]

Central first derivative, fourth order

(23)\[\left| \frac{\partial u}{\partial x}\right|_i = \frac{−u_{i+2}+8u_{i+1}−8u_{i−1}+u_{i−2}}{12\Delta x} +O(\Delta x)^4\]

Central second derivative, second order

(24)\[\left| \frac{\partial^2 u}{\partial x^2}\right|_i = \frac{u_{i+1}−2u_i+u_{i−1}}{\Delta x^2} +O(\Delta x)^2\]

Central second derivative, fourth order

(25)\[\left|\frac{\partial^2 u}{\partial x^2} \right|_i = \frac{−u_{i+2}+16u_{i+1}−30u_i+16u_{i−1}−u_{i−2}}{12\Delta x^2} +O(\Delta x)^4\]

Central third derivative, second order

(26)\[\left|\frac{\partial^3 u}{\partial x^3}\right|_i = \frac{u_{i+2}−2u_{i+1}+2u_{i−1}−u_{i−2}}{2\Delta x^3} +O(\Delta x)^2\]

Central third derivative, fourth order

(27)\[\left| \frac{\partial^3 u}{\partial x^3}\right|_i = \frac{−u_{i+3}+8u_{i+2}−13u_{i+1}+13u_{i−1}−8u_{i−2}+u_{i−3}}{8\Delta x^3} +O(\Delta x)^4\]

Central fourth derivative

(28)\[\left| \frac{\partial^4 u}{\partial x^4}\right|_i =\frac{u_{i+2}−4u_{i+1}+6u_i−4u_{i−1}+u_{i−2}}{\Delta x^4} +O(\Delta x)2\]

Note that the higher the order of the finite difference scheme, the more adjacent points are required. It is also important to note that derivatives of the following form

(29)\[\frac{\partial}{\partial x}\left(k\frac{\partial u}{\partial x}\right)\]

should be formed as follows

(30)\[\left| \frac{\partial}{\partial x}\left(k\frac{\partial u}{\partial x} \right)\right|_i = \frac{k_{i+1/2} \frac{u_{i+1}−u_i}{\Delta x} − k_{i−1/2} \frac{u_i−u_{i−1}}{\Delta x}}{\Delta x} +O(\Delta x)^2\]