Finite difference method, part II

All previous applications of finite difference methods to our problems have relied on the fact that we could reorganize the discretized equations in to form where the left

\[T_n^{i} = \textrm{some expression of } T^{i-1}\]

This is because we have chosen to discretize the equation with forward difference in time. For example with the simple heat diffusion equation (no advection):

\[\rho C_p \frac{T_n^{i} - T_n^{i-1}}{\Delta t} = \alpha \frac{T_{n+1}^{i-1} - 2T_n^{i-1} + T_{n-1}^{i-1}}{\Delta x^2} + H\]

We can, however, formulate the same heat equation with a backward difference in time (which is often beneficial for the stability of the solution):

(1)\[\rho C_p \frac{T_n^{i} - T_n^{i-1}}{\Delta t} = \alpha \frac{T_{n+1}^i - 2T_n^i + T_{n-1}^i}{\Delta x^2} + H\]

Now, in order to calculate the new value \(T_n^i\) we need to know the new values next to it, \(T_{n\pm 1}^i\). Clearly, this we are not able to use the new values to calculate those exact new values!

The way around this problem is to solve the system of linear equations formed by the difference equations. Every time we want to calculate the new values of \(T\) for the next time step, we have \(N\) equations, one for each grid point (\(x_0\), \(x_1\), \(x_2\), ..., \(x_{N-1}\)) Since this is also the number of unknown values of \(T\) s (\(T_0\), \(T_1\), \(T_2\), ..., \(T_{N-1}\)), we can solve all the values of \(T\) on one go.

Consider eq. (1). If we reorganize the equation to find all the coefficients for each \(T\), we get

\[T_{n-1}^i \left(\frac{-\alpha}{\Delta x^2}\right) + T_n^i \left(\frac{\rho C_p}{\Delta t} + \frac{2\alpha}{\Delta x^2}\right) + T_{n+1}^i \left(\frac{-\alpha}{\Delta x^2}\right) = T_n^{i-1} \rho C_p / \Delta t + H\]

Here, the coefficients for all the \(T\) s on the left hand side are known, and so is the whole right hand side, too. If we mark, for convenience, \(A=-\alpha/\Delta x^2\) and \(B=\frac{\rho C_p}{\Delta t} + \frac{2\alpha}{\Delta x^2}\), we can write

\[T_{n-1}^i A + T_n^i B + T_{n+1}^i A = T_n^{i-1} \rho C_p / \Delta t + H\]

For each spatial grid point \(n\) there exists one such equation. If we have four grid points in \(x\) direction, this leads a system of equations like

\[\begin{split}\begin{eqnarray} T_{-1}^i A +& T_0^i B +& T_1^i A & & & & = & T_0^{i-1} \rho C_p / \Delta t + H \\ & T_0^i A +& T_1^i B +& T_2^i A +& & & = & T_1^{i-1} \rho C_p / \Delta t + H \\ & & T_1^i A +& T_2^i B +& T_3^i A & & = & T_2^{i-1} \rho C_p / \Delta t + H \\ & & & T_2^i A +& T_3^i B +& T_4^i A & = & T_3^{i-1} \rho C_p / \Delta t + H \end{eqnarray}\end{split}\]

Since our grid has four points, there are no points \(n=4\) or \(n=-1\). So the first and last equation are not valid. We have to replace them by something given by the boundary conditions. If we have fixed surface temperature and fixed bottom temperature, we get following system:

\[\begin{split}\begin{eqnarray} T_0^i & & & & = & T_{\mathrm{surf}} \\ T_0^i A +& T_1^i B +& T_2^i A +& & = & T_1^{i-1} \rho C_p / \Delta t + H \\ & T_1^i A +& T_2^i B +& T_3^i A & = & T_2^{i-1} \rho C_p / \Delta t + H \\ & & & T_3^i & = & T_{\mathrm{bott}} \end{eqnarray}\end{split}\]

In matrix notation, we can state the same as

\[\begin{split}\begin{bmatrix} 1 & & & \\ A & B & A & \\ & A & B & A \\ & & & 1 \end{bmatrix} \begin{bmatrix} T_0^i \\ T_1^i \\ T_2^i \\ T_3^i \end{bmatrix} = \begin{bmatrix} T_{\mathrm{surf}} \\ T_1^{i-1} \rho C_p / \Delta t + H \\ T_2^{i-1} \rho C_p / \Delta t + H \\ T_{\mathrm{bott}} \end{bmatrix}\end{split}\]

or, in shorter notation

\[\mathbf{M}\mathbf{T}^i = \mathbf{b}\]

where \(\mathbf{M}\) is the coefficient matrix (values of which are known), \(\mathbf{T}^i\) is the temperature vector of unknown values (i.e. those we want to solve), and \(\mathbf{b}\) is the right hand side with known values.

In order to solve the unknown \(\mathbf{T}^i\) we can multiply the equation by the inverse of the matrix \(\mathbf{M}\).

\[\begin{split}\begin{eqnarray} \mathbf{M}^{-1} \mathbf{M} \mathbf{T}^i & = & \mathbf{M}^{-1} \mathbf{b} \\ \mathbf{T}^i & = & \mathbf{M}^{-1} \mathbf{b} \end{eqnarray}\end{split}\]

In Python (NumPy), there is a more direct way to achieve this, a function called np.linalg.solve().

See script heat_diff_simple_implicit.py for an example how to use it.

Exercise

Take the script template heat_diff_const_implicit.py and fill in the missing lines. Use heat_diff_simple_implicit.py as an example, if needed. This code will then calculate the cooling lithosphere problem using an implicit finite difference method.

Experiment with the size of time step (or, number of time steps taken): Can you make it “explode” like with the heat diffusion code with explicit formulation used previously?