Why does scipy Conjugate Gradient solver fail to converge for non-steady heat equation using Crank-Nicolson method

Could someone please explain why my implementation of the Crank-Nicolson method applied to the non-steady heat equation won't converge? There shouldn't be any nonlinear aspects to my implementation but scipy's sparse conjugate gradient solver won't seem to converge even when I ramp up the maxiter parameter to 1000. Scipy's least squares method finds a solution every time, seemingly with ease. Still, I would prefer to use scipy's indirect linear solvers if possible as I understand these to be more efficient for this type of problem.

I have included my implementation below. Any assistance is greatly appreciated. Thanks in advance!

EDIT The heat equation I am trying to solve is the following

\begin{align} \rho c \frac{\partial T}{\partial t} &= \frac{\partial}{\partial x}\big( \kappa (x) \frac{\partial T}{\partial x}\big) + q_{gen} (t) \; for \; (x,t) \in (0,1) \times (0,t_{dur}] \nonumber \\ T(x,0) &= T_0 \nonumber \\ \kappa (x) \frac{\partial T(0,t)}{\partial x} &= -\eta (T_{Amb}-T(0,t)) \nonumber \\ \kappa (x) \frac{\partial T(1,t)}{\partial x} &= \eta (T_{Amb}-T(1,t)) \end{align}

For thermal conductivity $$\kappa(x)$$, specific heat capacity $$c$$, density $$\rho$$, heat source $$q_{gen} (t)$$, and ambient temperature $$T_{Amb}$$. Note that while $$\kappa (x)$$ is constant in space in my below minimal-ish implementation, it is variable in my complete implementation.

Note also that I am using the Finite Volume (FV) method, hence the formulation in terms of flux $$q(x,t)$$ rather than temperature $$T(x,t)$$. Using notation from Patanka, S - Numerical heat transfer and fluid flow, the FV discretisation is found by integrating both sides of the PDE w.r.t. $$x$$ as per the following:

\begin{align} \int_w^e \rho \: c \frac{\partial T}{\partial t} \: dx &= \int_w^e \frac{\partial}{\partial x} \bigg( \kappa \frac{\partial T}{\partial x} \bigg) \: dx + \int_w^e q_{gen} \: dx \nonumber \nonumber \\ \rho c \frac{\partial T_P}{\partial t} \Delta x &= \bigg( \kappa \frac{\partial T}{\partial x } \bigg)_{e} - \bigg( \kappa \frac{\partial T}{\partial x } \bigg)_{w} + q_{gen} \Delta x \nonumber \end{align} Assuming temperature varies linearly between control volume faces \begin{align} \rho c \frac{\partial T_P}{\partial t} \Delta x &= \kappa_e \frac{(T_{E} - T_{P})}{\delta x_e} - \kappa_w \frac{(T_{P} - T_{W})}{\delta x_w} + q_{gen} \Delta x \nonumber \\ &= q_e - q_w + q_{gen} \Delta x \nonumber \\ \frac{\partial T_P}{\partial t} &= \frac{(q_e - q_w)}{\rho c \Delta x} + \frac{q_{gen}}{\rho c} \end{align}

Temperature is evaluated at grid point, $$P$$, with neighbours $$E$$ and $$W$$ (east and west respectively). Heat flux $$q(x)=\kappa (x) \frac{\partial T}{\partial x}$$ is evaluated at control volume faces denoted by $$e$$ and $$w$$ respectively. There are two spatial increments, control volume width $$\Delta x$$, and grid point spacing $$\delta x$$. These are both equal in my below implementation but not necessarily in my complete implementation. Half control volumes are used at the boundaries so the boundary flux can be used directly in the implementation.

As mentioned above, I am using the Crank-Nicolson (CN) method to discretise in time. Earlier iterations of my implementation adopted the explicit Euler method here but, to maintain stability, this was found to require time increments $$\Delta t$$ that were too small, thus making the total compute time impractical. The CN discretisation is as follows:

\begin{align} \frac{\partial T_P}{\partial t} &= \frac{(q_e - q_w)}{\rho c \Delta x} + \frac{q_{gen}}{\rho c} \nonumber \\ \frac{\partial T_P}{\partial t} &= F(x, t, q) \end{align}

For $$t^n = n \Delta t$$, $$F^n = F(x, t^n, q^n)$$ and $$F^{n+1} = F(x, t^{n+1}, q^{n+1})$$.

\begin{align} \frac{T^{n+1}_P - T^{n}_P}{\Delta t} &= \frac{F^{n+1} + F^{n}}{2} \nonumber \\ \frac{T^{n+1}}{\Delta t} - \frac{F^{n+1}}{2} &= \frac{T^{n}}{\Delta t} + \frac{F^{n}}{2} \nonumber \\ A(\vec{T^{n+1}}) &= \vec{b} \end{align}

from scipy.sparse import linalg
from scipy.optimize import least_squares
import numpy as np

# Mesh & time parameters
lx = 1.0
Nx = 10
tdur = 10
dt = 0.25

# Material properties
kappa = 9.7
rhoCp = 2200.
eta = 92.

# Heat source parameters
B = 0.012
C = 1.6
D = 6.2
t2 = 3.5
Hult = 380.
Cc = 500.

numGPs = Nx + 1
dx = lx/Nx
Np = int(tdur/dt)
tSeq = np.linspace(0, tdur, Np+1)
T0 = np.ones(numGPs)*25.
q0 = np.zeros(numGPs+1)
numSteps = int(tdur/dt)

# Ambient temperature parameters
tstart = 11     # Start time (h)
Tmin = 20.      # Min ambient temperature (C)
Tmean = 25.     # Mean ambient temperature (C)
Tmax = 30.      # Max ambient temperature (C)
TAmb = 0.5*(Tmin+Tmax) + np.sin(((tSeq + tstart - 9.)*360./24)*np.pi/180)*0.5*(Tmax - Tmin)

def logistic(L, x0, tSeq, k=500):
# Used as a "switch" for heat source equation
return L/(1 + np.exp(-k*(tSeq - x0)))

def get_H(B, C, D, t2, Hult, tSeq):
factor = logistic(0.5, t2, tSeq, 500)
H = 0.5*Hult*(1.0 - np.exp(-B*tSeq**C)) + factor*Hult*(tSeq - t2)/(tSeq - t2 + D)
return H

def get_qgen(B, C, D, t2, Hult, tSeq, Cc, dt):
H = get_H(B, C, D, t2, Hult, tSeq)
qgen = Cc*(H[1:]-H[:-1])/dt
return qgen

def get_F(Tn, kappa, rhoCp, eta, dx, TAmbn, qgenn):
dTdx = (Tn[1:] - Tn[:-1])/dx
qX = q0
qX[0] = -eta*(TAmbn - Tn[0])
qX[1:-1] = kappa*dTdx
qX[-1] = eta*(TAmbn - Tn[-1])
return ((qX[1:] - qX[:-1])/dx + qgenn)/rhoCp

def get_Tn(Tn_, Fn_, dt, kappa, rhoCp, eta, dx, TAmbn, qgenn):
"""
Tn_:    Temperature vector at time step, t=n.dt
Tn:     Temperature vector at time step, t=(n+1).dt
"""
def get_f_Tn(Tn):
Fn = get_F(Tn, kappa, rhoCp, eta, dx, TAmbn, qgenn)
return Tn/dt - Fn/2
b = Tn_/dt + Fn_/2

### Indirect linear solvers
A = linalg.LinearOperator((Tn_.shape[0], Tn_.shape[0]), matvec=get_f_Tn)
Tn, info = linalg.cg(A, b, Tn_)
#Tn, info = linalg.gmres(A, b, Tn_)
#Tn, info = linalg.bicgstab(A, b, Tn_)

### Nonlinear solvers
#result = least_squares(lambda Tn: get_f_Tn(Tn) - b, Tn_)
#Tn = result.x
#res = result.cost
#info = 0

res = np.linalg.norm(get_f_Tn(Tn) - b)
Fn = get_F(Tn, kappa, rhoCp, eta, dx, TAmbn, qgenn)
return Tn, Fn, res, info

Tn_ = T0
qgen = get_qgen(B, C, D, t2, Hult, tSeq, Cc, dt)
Fn_ = get_F(Tn_, kappa, rhoCp, eta, dx, TAmb[0], 0.)
for n, (TAmbn, qgenn) in enumerate(zip(TAmb[1:], qgen)):
Tn, Fn, res, info = get_Tn(Tn_, Fn_, dt, kappa, rhoCp, eta, dx, TAmbn, qgenn)
print('n=%3d, res=%.2e, info=%d'%(n+1, res, info))
#print(np.array2string(Tn, formatter={'float_kind':lambda x: "%.1f" % x}), '\n')
Tn_ = Tn
Fn_ = Fn
$$$$

• Is the matrix A symmetric or hermitian? This is often not the case if boundary conditions contain derivatives. Use gmres instead as solver. Commented Jul 22, 2023 at 12:57
• Also, basically all Krylov methods are guaranteed to converge in $N$ iterations, where $N$ is the size of the square input matrix. This is assuming the solver is being used on the correct class of matrices Commented Jul 22, 2023 at 13:24
• Could you write mathematically what you are solving, I cannot decipher your code. Are you solving $\partial_t u(t,x) = \Delta u(t,x)$ or $\partial_t u(t,x) = \operatorname{div}(\kappa(x) \nabla u(t,x)$? The latter would be unsymmetric, in which case you should use a solver like e.g. BiCGStab. If this is not the case have you tried implicit Euler? Commented Jul 22, 2023 at 14:54
• @lightxbulb ditto on the need for more exposition but you can see in get_F that there are some Robin conditions being applied consistent with Newton's law of cooling that will make the problem asymmetric. Now that there are some equations we can see that more clearly Commented Jul 22, 2023 at 22:41
• Thank you for your comments @LutzLehmann, @whpowell96 and @lightxbulb. I gave gmres and bicgstab a try without success - see updated implementation. It does sound like my implementation is asymmetric due to the boundary conditions though. Still, bicgstab` would be expected to converge? Is there a way to confirm the matrix class without actually instantiating the matrix? Also, I have tried implicit Euler but this becomes unstable for the range of material properties I am using. Commented Jul 22, 2023 at 22:49

Continuous Formulation

If I understood everything correctly you are trying to solve the heat equation on some domain $$\Omega\subseteq \mathbb{R}^d$$ with space-variant diffusivity $$\kappa:\Omega\to [0,\infty)$$, some term $$\rho:\Omega\to[0,\infty]$$, a forcing term $$q:\Omega\to\mathbb{R}$$, and Robin boundary conditions on the boundary $$\partial\Omega$$: \begin{alignat}{3} \newcommand{\bm}[1]{\boldsymbol{#1}} \rho(\bm{x})\partial_t T(t,\bm{x}) &= \bm{\nabla}\cdot (\kappa(\bm{x}) \bm{\nabla} T(t,\bm{x})) + q(\bm{x}), &\quad& \bm{x} \in \Omega, \,\,t\in(0,\infty), \\ \eta T(t,\bm{x}) + \kappa(\bm{x})\partial_{\bm{n}}T(t,\bm{x}) &= \eta T_{amb}, &\quad& \bm{x}\in\partial\Omega, \,\, t\in(0,\infty). \\ T(0,\bm{x}) &= T_0(\bm{x}), &\quad& \bm{x}\in\Omega. \end{alignat}

In the above $$\bm{n}:\partial\Omega\to\mathcal{S}^{d-1}$$ is the outward facing unit normal. For your specific case $$d=1$$, $$\Omega = (0,1)$$, $$\partial\Omega = \{0,1\}$$, $$n(0) = -1$$, $$n(1)=1$$, and thus $$\partial_{n}T(t,0) = -\partial_x T(t,0)$$ and $$\partial_n T(t,1) = \partial_x T(t,1)$$.

Weak Formulation

You want to discretise the problem with the finite volume method, so let us multiply both sides of the equation by a test function $$w:\Omega\to\mathbb{R}$$ and integrate over $$\Omega$$:

\begin{align} \int_{\Omega} \rho\, w\,\partial_t T &= \int_{\Omega} w\,\bm{\nabla}\cdot (\kappa\bm{\nabla} T) + \int_{\Omega}w\,q \\ &= \int_{\Omega}\bm{\nabla}\cdot(w\kappa\bm{\nabla}T) - \int_{\Omega}\kappa\,\bm{\nabla}w \cdot \bm{\nabla}T + \int_{\Omega} w \,q \\ &=\int_{\partial\Omega}w\kappa\,\partial_{\bm{n}}T -\int_{\Omega}\kappa\,\bm{\nabla}w \cdot \bm{\nabla}T + \int_{\Omega} w \,q. \end{align}

FVM discretisation: Piecewise Constant Test Functions

Now consider a finite-dimensional space for the test functions made up of piecewise-constant functions on a partition of $$\Omega$$ made up of the control volumes $$\{K_0,\ldots,K_{m+1}\}$$:

$$$$W = \left\{w(x) = \sum_{j=0}^{m+1} \alpha_j \bm{1}_{K_j}(x)\right\}, \, \Omega = \bigcup_{i=0}^{m+1} K_i, \, K_i\ne \emptyset, \, K_i \cap K_j = \emptyset, \, i\ne j.$$$$

A basis of that space are of course the indicator functions: $$\psi_i(\bm{x}) = \bm{1}_{K_i}(\bm{x})$$. Then instead of integrating over $$\Omega$$ we may integrate over $$K_i$$ and set $$w=\psi_i$$ which yields:

$$$$\int_{K_i}\rho \partial_tT = \int_{\partial K_i}\kappa\, \partial_{\bm{n}} T + \int_{K_i} q,$$$$

where the term $$\int_{K_i} \kappa \, \bm{\nabla} w \cdot \bm{\nabla} T = 0$$ since $$\bm{\nabla} \psi_i = 0$$.

Partition in 1-D

In your 1-D setting let the faces be given as $$0=f_0, i.e.\ $$K_i = (f_i,f_{i+1})$$ and let the centers $$(x_i)_{i=0}^{m+1}$$ be such that $$x_i\in \overline{K_i}$$, where $$\overline{K_i} = [f_i,f_{i+1}]$$ (note that you should probably avoid degenerate cases where $$x_i=x_{i+1}$$). If you don't have any specific knot vector for the $$x_i$$ you may set them at the centers of the control volumes: $$x_{i+1} = \frac{f_i+f_{i+1}}{2}$$. For the discretisation I will use, the boundary points should be chosen such that $$0\leq i \leq m-1$$, and $$x_0 = f_0 = 0$$, $$x_{m+1} = f_{m+2}$$ = 1.

Flux Discretisation: Continuous Piecewise Linear Trial Functions

We would like to approximate the flux terms $$\kappa(f_i)\partial_xT(f_i)$$. I have chosen to approximate the solution $$T$$ with linear function in every interval $$[x_i,x_{i+1}]$$ and require that it is continuous on $$\Omega$$. Note that this choice is consistent with the discretisation you get from finite differences on an irregular grid, the one you get from the finite element method with linear elements, and the one you get from the finite volume method with the two-point flux approximation scheme. Our solution resides in the space of continuous piecewise linear functions (i.e.\ polygonal chains) on the knot vector $$(x_0,\ldots,x_{m+1})$$:

$$$$V = \left\{v(\bm{x}) = \sum_{j=0}^{m+1}\alpha_j \phi_j(\bm{x})\right\},$$$$

where the basis functions $$\phi_i$$ are the hat functions:

$$$$\phi_i(x) = \begin{cases} \frac{x-x_{i-1}}{x_i-x_{i-1}}, &\text{for } x \in K_{i-1}, \\ \frac{x-x_{i+1}}{x_{i}-x_{i+1}}, &\text{for } x \in K_i, \\ 0, &\text{for } x\in\Omega \setminus (K_{i-1} \cup K_i). \end{cases}$$$$

Since the functions are linear and $$x_i\leq f_i \leq x_{i+1}$$ we have: $$\partial_{x}T(f_i) = \frac{T(x_{i})-T(x_{i-1})}{x_{i}-x_{i-1}}=(T_{i} - T_{i-1}) \Delta_{i-1}$$ for $$1\leq i \leq m$$.

You mention that $$\kappa$$ is discontinuous at faces. One option would be to take:

$$$$\kappa_{i-1/2} = \frac{\lim_{h\to 0^+}\kappa(f_i+h) + \lim_{h\to 0^+}\kappa(f_i-h)}{2}.$$$$

If $$\kappa$$ is a constant function on the volumes then that would be just the average between $$\kappa$$ on $$K_{i-1}$$ and on $$K_{i}$$. If you want a more formal treatment of this then you should probably look into discontinuous Galerkin methods. If you assume instead that $$\kappa$$ varies linearly between the control volume centers, you can take the linearly interpolated value: $$$$\kappa_{i-1/2} :=\kappa(f_i) \approx \frac{f_i-x_i}{x_{i-1}-x_i}\kappa(x_i) + \frac{f_i-x_{i-1}}{x_{i}-x_{i-1}}\kappa(x_{i-1}).$$$$

I don't know which you should take - you could try both.

Control Volume Integral Discretisation

The simplest discretisation of the control volume integrals is by assuming that the functions $$\rho, q, T$$ are constant on the cells $$K_i$$, then:

$$$$\int_{K_i}\rho\partial_t T \approx |K_i| \rho(x_i) \partial_t T(x_i), \quad \int_{K_i} q \approx |K_i| q(x_i), \quad |K_i| = f_{i+1}-f_i.$$$$

As a "better" alternative you may replace $$\rho(x_i)$$ with $$\rho_{avg,i} = \int_{K_i}\rho / |K_i|$$, and similarly for $$q(x_i)$$.

You may rightfully object that just before that we claimed that $$T$$ belongs to a space of continuous piecewise linear functions. In fact this mix of approximations is even called the finite volume element method sometimes for that reason (also box method). In finite element literature they call the resulting matrix a lumped mass matrix, it avoids having to solve an additional linear system since the mass matrix is diagonal in that case.

Of course, if you wish you may discretise $$\partial_t T$$ as a linear function on $$K_i$$ and integrate, then the result would also depend on the the two neighbours neighbours.

Space-discrete Formulation

Substituting our discretisation in the weak formulation projected onto the space of piecewise constant test functions yields:

\begin{align} \int_{K_i}\rho \partial_tT &= \int_{\partial K_i}\kappa\, \partial_{\bm{n}} T + \int_{K_i} q, \\ \rho_{avg,i} |K_i|d_tT_i &= \kappa_{i+1/2}(T_{i+1}-T_i)\Delta_i - \kappa_{i-1/2}(T_i-T_{i-1})\Delta_{i-1} + |K_i| q_{avg,i}, \,\, 1\leq i \leq m. \end{align}

Note that this produced $$m$$ equations for the volumes $$K_1,\ldots,K_m$$. We also need to handle the boundary conditions and boundary volumes: \begin{align} \eta T(t,\bm{x}) + \kappa(\bm{x})\partial_{\bm{n}}T(t,\bm{x}) &= \eta T_{amb}, \quad \bm{x}\in\partial\Omega, \,\, t\in(0,\infty), \\ \eta_0 T_0 - \kappa_0 (T_1-T_0)\Delta_0 &= \eta_0 T_{amb}, \\ \eta_{m+1} T_{m+1} + \kappa_{m+1} (T_{m+1}-T_m)\Delta_m &= \eta_{m+1} T_{amb}, \end{align}

You can use the above to solve for $$T_0$$ and $$T_{m+1}$$ as:

\begin{align} T_0 = \frac{\eta_0 T_{amb} + \kappa_{0}\Delta_0 T_1}{\eta_0+\kappa_0\Delta_0}, \quad T_{m+1} = \frac{\eta_{m+1}T_{amb} - \kappa_{m+1}\Delta_mT_{m+1}}{\eta_{m+1}-\kappa_{m+1}\Delta_m}. \end{align}

You may then substitute this solution for $$T_0$$ in the equation for $$d_tT_1$$, and $$T_{m+1}$$ in the equation for $$d_tT_m$$, which would eliminate the variables and you would be able to work with an $$m\times m$$ matrix. Then at the end you can compute $$T_0$$ and $$T_{m+1}$$ from the solution for $$T_1$$ and $$T_m$$. The other option is of course to keep these two rows, then you would have an $$(m+2)\times (m+2)$$ matrix. In either case we may write the system in matrix-vector form as:

$$$$\bm{M}_{\bm{\rho}}d_t\bm{T}(t) = -\bm{W}_{\bm{\kappa}}\bm{T}(t) + \bm{q},\quad \bm{T}(0) = \bm{T}_0.$$$$

Time-continuous Solution

If we assume that $$\rho_{avg,i}$$ and $$|K_i|$$ are non-zero for all $$i$$, then the matrix $$\bm{M}_{\bm{\rho}}$$ is invertible, and we may write:

\begin{align} d_t\bm{T} &= -\bm{M}_{\bm{\rho}}^{-1}\bm{W}_{\bm{\kappa}}\bm{T}(t) + \bm{M}_{\bm{\rho}}^{-1}\bm{q} = -\bm{L}\bm{T}(t) + \bm{b},\\ \bm{T}(t) &= \bm{\exp}(-t\bm{L})\bm{T}_0 + \left(\int_{0}^t \bm{\exp}((s-t)\bm{L})\,ds\right)\bm{b}, \\ &= \bm{\exp}(-t\bm{L})\bm{T}_0 + \left(\sum_{k=0}^{\infty} \frac{t(-t\bm{L})^k}{(k+1)!}\right)\bm{b}, \\ &= \bm{\exp}(-t\bm{L})\bm{T}_0 + \biggl(\bm{L}^{\#} - \bm{\exp}(-t\bm{L})\bm{L}^{\#}+t(\bm{I}-\bm{L}\bm{L}^{\#})\biggr)\bm{b}, \end{align}

where $$\bm{L}^{\#}$$ is the Drazin inverse. This solution doesn't require a discretisation in time, but instead requires computing matrix exponentials, and this Drazin inverse. While both are easy to do if you have the eigendecomposition $$\bm{L} = \bm{P}\bm{\Lambda}\bm{P}^{-1}$$, it is not so easy otherwise. However both explicit Euler and Implicit Euler approximate this solution as the step size goes to zero.

Series solution

One option is to directly use one of the solutions from earlier and write the exponential as a series:

\begin{align} \bm{T}(t) &= \bm{\exp}(-t\bm{L})\bm{T}_0 + \left(\sum_{k=0}^{\infty} \frac{t(-t\bm{L})^k}{(k+1)!}\right)\bm{b} \\ &= \left(\sum_{k=0}^{\infty}\frac{(-t\bm{L})^k}{k!}\right) \bm{T}_0 + \left(\sum_{k=0}^{\infty} \frac{t}{k+1}\frac{(-t\bm{L})^k}{k!}\right)\bm{b} \\ &= \sum_{k=0}^{\infty}\frac{(-tL)^k}{k!} \left(\bm{T}_0 + \frac{t}{k+1}\bm{b}\right). \end{align}

I can't tell you when you can truncate this series, but it should give you an approximation of the solution provided you truncate after enough terms. You don't have to compute powers of matrices either, you can just perform two matrix-vector mutliplications to update some product variables. If you can compute the Drazin inverse applied to $$\bm{b}$$ (e.g.\ using some iterative method) then you can reduce that to one matrix-vector multiplication :

\begin{align} \bm{T}(t) &= \bm{\exp}(-t\bm{L})\bm{T}_0 + \biggl(\bm{L}^{\#} - \bm{\exp}(-t\bm{L})\bm{L}^{\#}+t(\bm{I}-\bm{L}\bm{L}^{\#})\biggr)\bm{b} \\ &= \sum_{k=0}^{\infty}\frac{(-tL)^k}{k!}(\bm{T}_0-(\bm{L}^{\#}\bm{b})) + (\bm{I}+t\bm{I}-t\bm{L})(\bm{L}^{\#}\bm{b}). \end{align}

I believe that Matlab uses a Pad{'e} approximation of the exponential for faster series convergence.

Explicit Euler

Explicit Euler can be formulated as taking forward differences in time for $$d_t\bm{T}$$, and taking the right-hand side at the old time step:

$$$$\frac{\bm{T}^{k+1}-\bm{T}^k}{\tau} = -\bm{L}\bm{T}^k + \bm{b} \implies \bm{T}^{k+1} = (\bm{I}-\tau\bm{L})\bm{T}^k + \tau\bm{b}.$$$$

Let the steady state solution be $$\bm{T}^*$$ (this is the solution at time infinty), then the error is given as

\begin{align} \bm{\epsilon}^{k+1} &= \bm{T}^{k+1}-\bm{T}^* = (\bm{I}-\tau\bm{L})\bm{T}^k + \tau \bm{b} - (\tau\bm{b}-\tau\bm{L}\bm{T}^*) \\ &= (\bm{I}-\tau\bm{L})(\bm{T}^k-\bm{T}^*) = (\bm{I}-\tau\bm{L})\bm{\epsilon}^k \\ &= (\bm{I}-\tau\bm{L})^{k+1}\bm{\epsilon}^0 = (\bm{I}-\tau\bm{L})^{k+1}(\bm{T}_0-\bm{T}^*). \end{align}

Since we want for the error to go to zero $$\lim_{k\to\infty}\|\bm{e}^k\| = \lim_{k\to\infty}\|(\bm{I}-\tau\bm{L})^k\bm{e}^0\| = 0$$ for some norm $$\|\cdot\|$$, and knowing that $$\|\bm{e}^{k+1}\| = \|(\bm{I}-\tau\bm{L})\bm{e}^k\|\leq \|\bm{I}-\tau\bm{L}\|\|\bm{e}^k\|$$, where $$\|\cdot\|$$ is the matrix norm induced by $$\|\cdot\|$$, a sufficient condition is $$\|\bm{I}-\tau\bm{L}\|<1$$. For the $$2$$-norm and a symmetric positive semi-definite matrix $$\bm{L}$$ the condition is then $$|1-\tau\lambda_i|<1$$ for any $$i$$, which yields $$\tau\in (0, 2/\lambda_{\max})$$. This is typically where your step size constraint comes from.

However in your case you have a non-symmetric matrix, so the $$2$$-norm is given by the largest singular value of $$(\bm{I}-\tau\bm{L})$$, thus the constraint reads:

$$$$\|(\bm{I}-\tau\bm{L})^T(\bm{I}-\tau\bm{L})\| = \|\bm{I}-\tau(\bm{L}^T+\bm{L}) + \tau^2 \bm{L}^T\bm{L}\|<1.$$$$

So I am not even sure you can guarantee there always exists a step size that works for non-symmetric matrices. It is for instance clear that this method doesn't work even for symmetric indefinite matrices, because you can never find a step size that would decrease the error for all possible $$\bm{e}^k$$, since for negative eigenvalues you get $$1-\tau \lambda_i < 1 \implies \tau<0$$, while for positive eigenvalues you have $$\tau>0$$, i.e. there is no such $$\tau$$. Note that explicit Euler is equivalent to the Richardson iteration (also to Jacobi if you divide by the diagonal).

Implicit Euler

Implicit Euler discretises the time derivative with a forward difference but takes the right hand side at the future time:

$$$$\frac{\bm{T}^{k+1}-\bm{T}^k}{\tau} = -\bm{L}\bm{T}^{k+1} + \bm{b} \implies (\bm{I}+\tau\bm{L})\bm{T}^{k+1} = \bm{T}^k+\tau\bm{b}.$$$$

Note that you have to solve a linear system for each step of the above, in your case this system is non-symmetric so you should use BiCGStab, GMREs, CGNE, CGNR/LSQR for it. Since your system is tridiagonal you can also use the Thomas algorithm directly. If $$\bm{L}$$ is symmetric positive semi-definite then this is unconditionally stable, i.e. you can take any $$\tau > 0$$. For the indefinite symmetric case implicit Euler has the constraint $$\tau > \frac{2}{-\lambda_{neg,\min}}$$, where $$\lambda_{neg,\min}$$ is the negative eigenvalue that is closest to zero, so you can also construct examples where it is unstable. That is, implicit Euler can actually be unstable because you took a step that is too small if your matrix is symmetric and has negative eigenvalues. This makes sense because the iteration matrix is $$(\bm{I}+\tau\bm{L})^{-1}$$, thus the eigenvalues are reciprocal, and you get the condition $$|(1+\tau\lambda_i)^{-1}|<1$$. For non-symmetric matrices you need to again consider the 2-norm of the iteration matrix, and then you get

$$$$\|(\bm{I}+\tau\bm{L})^{-T}(\bm{I}+\tau\bm{L})^{-1}\| = \|(\bm{I}+\tau(\bm{L}^T+\bm{L}) + \tau^2\bm{L}^T\bm{L})^{-1}\|<1$$$$

Crank-Nicolson

The Crank-Nicolson discretisation has the form:

\begin{aligned} \bm{T}^{k+1}-\bm{T}^k &= -\frac{\tau}{2}\bm{L}\bm{T}^k - \frac{\tau}{2}\bm{L}\bm{T}^{k+1} + \tau\bm{b}, \\ (\bm{I}+0.5\tau\bm{L})\bm{T}^{k+1} &= (\bm{I}-0.5\tau\bm{L})\bm{T}^k + \tau\bm{b}, \\ \bm{T}^{k+1} &= (\bm{I}+0.5\tau\bm{L})^{-1}(\bm{I}-0.5\tau\bm{L})\bm{T}^k + \tau(\bm{I}+0.5\tau\bm{L})^{-1}\bm{b}. \end{aligned}

Once again you have to solve a linear system of equations to get the solution. The supposed benefit is that your error decays as $$O(\tau^2)$$ in time. For symmetric positive semi-definite matrices $$\bm{L}$$ this is again unconditionally stable for $$\tau>0$$ since we have:

$$$$\|(\bm{I}+0.5\tau\bm{L})^{-1}(\bm{I}-0.5\tau\bm{L})\|<1 \iff \left|\frac{1-0.5\tau \lambda_i}{1+0.5\tau\lambda_i}\right| < 1.$$$$

The latter holds for any $$\tau>0$$, but breaks if you have even one zero or negative eigenvalue. The zero case may be admitted however, as we may not want to remove the error corresponding to a zero eigenvector (e.g. if we wish to preserve the average grey value). Once again for unsymmetric matrices you have the condition:

$$$$\|(\bm{I}-0.5\tau\bm{L})^T(\bm{I}+0.5\tau\bm{L})^{-T}(\bm{I}-0.5\tau\bm{L})^{-1}(\bm{I}-0.5\tau\bm{L})\|<1,$$$$

which is again hard to argue about. Note that Crank-Nicolson benefits you if you want to take small time steps because of the $$O(\tau^2)$$ error decay. However if you plan to take large time steps, e.g. larger than the bound for explicit Euler, then you may get nasty manifestations of the error in the solution, from wikipedia:

However, the approximate solutions can still contain (decaying) spurious oscillations if the ratio of time step $$\Delta t$$ times the thermal diffusivity to the square of space step, $$(\Delta x)^2$$, is large (typically, larger than $$1/2$$ per Von Neumann stability analysis).

A Warning

Note that I haven't tested any of the above nor have I verified it extensively, so there may be typos and mistakes. Also I am not sure whether I understood correctly the formulation that you were going for.

• Brilliant! Thanks very much @lightxbulb. I'll go through your response properly when I finish work. One quick comment though, and apologies for not making this clear in my post, in my complete implementation I have material boundaries that align with control volume faces such that each control volume only has one set of material properties. This is the main driver for calculating the flux at the control volume faces and therefore assuming that temperature varies linearly between control volume faces. Commented Jul 24, 2023 at 3:01
• @n1ck94 What do you mean by material properties? The diffusivity $\kappa$? Or do you mean that $\rho$ and $c$ are also space variant? Commented Jul 24, 2023 at 7:31
• i'm aware that we are venturing into territory that goes beyond the scope of my initial post but in my complete implementation, each of $\kappa$, $\rho$ and $c$ vary in space and have discontinuities at material interfaces. Also, not all control volumes are the same width (due to each material having a different width), i.e. the face location $f_i$ is not always the average of the two surrounding centres. I'm happy to create separate posts to tackle the space-variant $\rho$ and $c$, and the variable mesh. Can discontinuities in $\kappa$ at control volume faces be accommodated? Commented Jul 24, 2023 at 12:24
• Also, it is very interesting that the explicit Euler method is possibly stable for any $\tau > 0$, this would be my preferred method for discretising time. I'm curious as to why the von Neumann numerical stability condition does not apply though (see the Stability section on Wikipedia page en.wikipedia.org/wiki/FTCS_scheme). Would you mind elaborating on this? Thanks! Commented Jul 24, 2023 at 12:37
• @n1ck94 I extended it to handle space varying $\rho$ can $c$, I combined both of them into a single function $\rho$ for conciseness though, since they are grouped anyways. Discontinuities in $\kappa$ at the faces are possible, and I think even what I suggested should work. I am not sure however, you can look into discontinuous Galerkin methods for more details - I linked a good book. It's not explicit Euler that is unconditionally stable, it is implicit Euler. And by unconditionally one really means for symmetric/Hermitian positive semi-definite matrices. I have sketched some of that out above Commented Jul 24, 2023 at 17:05