Implicit Diffusion and the Thomas Algorithm — Buying Stability with One Tridiagonal Solve

Halving the mesh quarters the explicit time step. In one dimension the way out is a single tridiagonal sweep, not a smaller step.

The run left on at 11pm was NaNs by morning. The grid had just been cut in half. The time step had not. Forgetting the CFL (Courant–Friedrichs–Lewy) line cost an hour of cluster time.

Halve the spacing and explicit diffusion quarters the step. A 1024-point mesh still walks hundreds of thousands of steps, only so it will not explode. The exit is short. Put the future on the right-hand side, then let Thomas factor the tridiagonal.

The r slider runs from 0.05 to 5 and lays the same Gaussian on both schemes, so the split shows up on one plot.

CFL number r = D·dt / dx²: 0.45Explicit: stable (r ≤ 0.5)Implicit: unconditionally stable

Below 0.5 the curves sit on top of each other. Cross 0.5 and the orange explicit trace starts to sawtooth; past 1 it leaves the frame, while the cyan implicit trace stays put.

CFL Locks Down Fine Grids#

The diffusion equation is

\frac{\partial q}{\partial t} = D \frac{\partial^2 q}{\partial x^2}

where $q$ is concentration (or temperature) and $D$ the diffusion coefficient. The simplest explicit scheme (forward Euler with central differences) reads

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n} - 2 q_i^{n} + q_{i-1}^{n}),\quad r = \frac{D\,\Delta t}{\Delta x^2}

For stability we need $r \le 1/2$ . A one-line von Neumann analysis: the shortest wavelength (two cells) has amplification factor $|1 - 4r|$ , and forcing it $\le 1$ gives $r \le 1/2$ .

Going from $1024 \to 2048$ cuts $\Delta t$ by four. One more cut, to $4096$ , makes $\Delta t$ one-sixteenth of the 1024-point step, so the step count grows by 16. Eight times the points means 64 times the steps to reach the same instant. The bill is that square, not the cell count.

At $r = 0.6$ the shortest-wave factor is $|1 - 4r| = 1.4$ . Sixty steps is that factor raised to the 60th power, eating the mode. The printed $1.33\times 10^{12}$ also carries whatever piece of the initial pulse sat on that mode. At $r = 0.4$ the factor is $0.6$ , and the explicit peak stays at $4.21\times 10^{-1}$ .

Push the Future to the Right Side#

One substitution. Replace $q^n$ on the right with $q^{n+1}$ .

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n+1} - 2 q_i^{n+1} + q_{i-1}^{n+1})

The future value sits on the right, so the cells cannot be updated one at a time. Every unknown is tied to its neighbors, and the step is a matrix equation.

A\,\mathbf{q}^{n+1} = \mathbf{q}^{n}

The same analysis gives an amplification factor $1/(1 + 4r)$ , always below 1 for positive $r$ . Raising $\Delta t$ does not make the scheme diverge. It is unconditionally stable.

Accuracy is the price. A huge $\Delta t$ stays stable and stops being accurate. "Will not blow up" is still worth a lot when the job runs overnight.

At the same $r = 0.4$ and the same 60 steps, the peaks are $4.21\times 10^{-1}$ explicit and $0.4456$ implicit. Both are stable, and they already disagree in the third digit. Same physical time. The gap is truncation error.

Tridiagonal Matrices and Thomas's Single Sweep#

Write $A$ . In 1D, cell $i$ couples only to $i-1$ and $i+1$ .

A = \begin{pmatrix} 1+r & -r & & & \\ -r & 1+2r & -r & & \\ & -r & 1+2r & -r & \\ & & \ddots & \ddots & \ddots \\ & & & -r & 1+r \end{pmatrix}

This is a tridiagonal matrix. Only three diagonals are nonzero. Generic Gaussian elimination is $\mathcal{O}(N^3)$ ; a tridiagonal system is $\mathcal{O}(N)$ . The algorithm is Thomas's, written up by Llewellyn Thomas as an IBM memo in 1949.

Thomas makes two passes.

Forward elimination — sweep top to bottom, eliminating the sub-diagonal while computing a modified diagonal $b'$ and right-hand side $d'$ .
Backward substitution — sweep bottom to top, plugging in known values.

Each pass is four or five flops per cell, and storage is three arrays. The full matrix never has to sit in memory.

Only the boundary diagonal is $1+r$ . Interior rows are $1+2r$ . Zero flux sets the ghost cell equal to the boundary value, so the outside neighbor drops out. One $r$ leaves the $2r$ on the diagonal, and the $-r$ that pointed at that neighbor leaves with it. That is why the code sets b[0] = 1 + r. Copy the interior row onto the boundary and the total leaks.

Comparing in NumPy#

import numpy as np
 
def explicit_diffusion_step(q, r):
    """Forward Euler + central difference."""
    next_q = q.copy()
    next_q[1:-1] = q[1:-1] + r * (q[2:] - 2 * q[1:-1] + q[:-2])
    return next_q
 
def thomas_solver(a, b, c, d):
    """Solve tridiagonal system A x = d in O(N).
    a: sub-diagonal (a[0] is unused)
    b: main diagonal
    c: super-diagonal (c[-1] is unused)
    d: right-hand side
    """
    n = len(d)
    cp = np.zeros(n)
    dp = np.zeros(n)
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / m
        dp[i] = (d[i] - a[i] * dp[i - 1]) / m
    x = np.zeros(n)
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def implicit_diffusion_step(q, r):
    """Backward Euler — solve (I - r·L) q^{n+1} = q^n via Thomas."""
    n = len(q)
    a = np.full(n, -r)
    b = np.full(n, 1 + 2 * r)
    c = np.full(n, -r)
    b[0] = 1 + r        # zero-flux boundary
    b[-1] = 1 + r
    return thomas_solver(a, b, c, q.copy())
 
# Compare on a Gaussian pulse
N = 80
x = np.linspace(0, 1, N)
q0 = np.exp(-((x - 0.5) ** 2) / 0.005)
D, dx = 1.0, 1.0 / N
r_values = [0.4, 0.6, 2.0]
for r in r_values:
    dt = r * dx**2 / D
    qe, qi = q0.copy(), q0.copy()
    for _ in range(60):
        qe = explicit_diffusion_step(qe, r)
        qi = implicit_diffusion_step(qi, r)
    print(f"r={r:.2f} | explicit max={np.max(np.abs(qe)):.2e} | implicit max={np.max(np.abs(qi)):.4f}")
# r=0.40 | explicit max=4.21e-01 | implicit max=0.4456
# r=0.60 | explicit max=1.33e+12 | implicit max=0.3873
# r=2.00 | explicit max=2.78e+38 | implicit max=0.2517

At $r=0.6$ the explicit method has exploded to $10^{12}$ in 60 steps. Implicit, on the same step, keeps spreading, with a peak of $0.3873$ . At $r=2$ the explicit max is $2.78\times 10^{38}$ and the implicit peak has fallen to $0.2517$ .

The three printed rows are not the same time#

With $\Delta t = r\,\Delta x^2 / D$ , $N = 80$ , and $D = 1$ , sixty steps of physical time scale with $r$ . The $r = 2$ row is five times farther along than the $r = 0.4$ row. Reading $0.4456$ next to $0.2517$ and calling the implicit scheme "more dead" is a false comparison. Matching the time means matching $r$ times the step count. Twelve steps at $r = 2$ land on the same instant as sixty steps at $r = 0.4$ . The printout never made that match.

Why Implicit Isn't Always the Answer#

Unconditional stability is not a free pass. The cost splits three ways.

1. Accuracy vs. stability — Stability says the amplitude will not explode. It does not say the answer is right. Backward Euler is $\mathcal{O}(\Delta t)$ . When the large step also has to be accurate, Crank–Nicolson ( $\mathcal{O}(\Delta t^2)$ , still unconditionally stable) is the better switch.

The shortest-wave factor for Crank–Nicolson is $(1-2r)/(1+2r)$ . At $r = 2$ that is $-0.6$ . The magnitude stays under 1, so it does not blow up, but the sign is negative and the phase flips every step. Push $r$ higher and the factor sticks to $-1$ .

2. Nonlinearity — If $D$ depends on the solution ( $D = D(q)$ ), the matrix changes every step and a nonlinear solver (Picard or Newton) has to sit on top of Thomas.

3. Higher dimensions — In 2D/3D the matrix is no longer tridiagonal. A five- or seven-point stencil grows side bands. Split the step with ADI (Alternating Direction Implicit) into two or three 1D tridiagonal sweeps, or move to iterative methods from Jacobi through BiCGSTAB.

That third case shows up again on the way from compressible to incompressible flow. The pressure Poisson equation $\nabla^2 P = -\rho \nabla\!\cdot(\mathbf{u}\!\cdot\!\nabla\mathbf{u})$ has the same large sparse shape. Anyone who has already solved implicit diffusion uses the same tools on the pressure step.

Calling a large r the right answer is a failure#

The implicit peak $0.2517$ at $r = 2$ is not a record of getting closer to the truth. Physical time is five times longer, and the time error scales with $\Delta t$ . Saving the curve without writing down the physical time of those steps makes the comparison a lie.

Once the grid is two-dimensional this tridiagonal code cannot be reused as it stands. When a sweep that never forms the matrix looks cheaper, how far does LU-SGS, which advances implicitly without building the matrix carry the same stability promise? If the sweep direction is wrong, does that promise break with it? That part is still open.