The Thomas Algorithm

The Thomas algorithm solves tridiagonal systems in linear time, the specialized elimination behind splines and one-dimensional implicit schemes.

A tailored elimination

Many problems produce tridiagonal systems, where each equation involves only an unknown and its two immediate neighbors. The Thomas algorithm is Gaussian elimination specialized to this structure: it performs a forward sweep to eliminate the subdiagonal and a backward sweep to substitute, all in O(n) operations rather than O(n^3).

def thomas(a, b, c, d):
    # a: sub, b: diag, c: super, d: rhs (a[0], c[-1] unused)
    n = len(b); cp = [0.0]*n; dp = [0.0]*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 = [0.0]*n; x[-1] = dp[-1]
    for i in range(n-2, -1, -1):
        x[i] = dp[i] - cp[i]*x[i+1]
    return x

Where it appears

Tridiagonal systems arise constantly: cubic spline coefficients, one-dimensional finite-difference discretizations of second derivatives, and the implicit steps of the Crank-Nicolson and ADI schemes for parabolic PDEs. In each case the Thomas algorithm gives an exact, fast, and stable solve.

Stability

The algorithm is stable without pivoting when the matrix is diagonally dominant or symmetric positive definite, which covers most physical discretizations. If those conditions fail, pivoting may be needed, slightly complicating the sweep. For the common well-conditioned cases, the plain algorithm is both fast and accurate.

Extensions

Block-tridiagonal systems, where each entry is itself a small matrix, generalize the method for coupled one-dimensional problems. Cyclic reduction and parallel variants break the sequential dependency for many-core machines. The banded generalization handles a few extra diagonals with the same linear scaling.

Tridiagonal solves appear in line-implicit and directional-splitting schemes used within breeder Hyperion transport solvers.