Tridiagonal matrix algorithm ($p_i=b_i-a_i{c}_{i-1}/p_{i-1}$)

The Thomas algorithm is Gaussian elimination specialized to a tridiagonal matrix: each row couples an unknown only to its immediate neighbours. Eliminating one unknown changes just the next diagonal entry and right-hand side. Provided all pivot elements are nonzero, a forward sweep followed by back substitution solves an $n$-unknown linear system in $O(n)$ time and storage, using Big O notation. Treating the same entries as a full matrix wastes the zeros: classical dense LU factorization and solve costs grow cubically for factorization and quadratically for storage.

A concrete reason to care is the Dirichlet Poisson solution by tridiagonal elimination: a three-point finite difference for the second derivative couples adjacent points on a mesh, producing precisely this tridiagonal matrix. The same structure occurs in one-dimensional Crank-Nicolson method calculations. NPTEL's archived Lecture 13 introduces the Thomas algorithm from that time-stepping problem and labels the procedure “Thomas algorithm (1949)”.^[ref]

Let $n\ge 1$. Store the tridiagonal matrix $A$ by its diagonal $b_1,\dots ,b_n$, lower diagonal $a_2,\dots ,a_n$, and upper diagonal $c_1,\dots ,c_{n-1}$. The linear system $Ax=d$ iswith the nonexistent endpoint terms omitted. Initialize $p_1=b_1$ and $g_1=d_1$. For $i=2,\dots ,n$, compute the elimination multiplier, modified pivot element, and modified right-hand side:These recurrence relations apply the elementary row operation $R_i\leftarrow R_i-m_i{R}_{i-1}$ after the preceding row has already been reduced. The transformed linear system is an upper bidiagonal matrix:Thus back substitution is especially short:For $n=1$, only $x_1=d_1/b_1$ remains. Check each pivot element before its first division, including $p_n$ before back substitution.

Equivalently, this constructs an LU decomposition with unit lower bidiagonal matrix $L$, whose lower diagonal is $m_2,\dots ,m_n$, and upper bidiagonal matrix $U$, whose diagonal is $p_1,\dots ,p_n$ and upper diagonal is $c_1,\dots ,c_{n-1}$. The $g$ sweep is forward substitution in $Lg=d$. There is no fill-in outside these bands when no rows are interchanged. LAPACK's DGTSV source uses these same coefficient and right-hand-side updates in its branch without row interchange.^[ref]

Consider the original small Poisson equation exampleIts exact solution is the quadratic polynomial $u(x)=1+2x-x^2$. Use a uniform mesh $x_i=ih$ with $h=1/4$ and three interior unknowns $U_1,U_2,U_3$. The finite difference method replaces the negative second derivative byMultiplying by $h^2$ and moving the prescribed Dirichlet boundary conditions to the right givesThe endpoint values contribute $1$ and $2$ to the first and last equations; they are not additional unknowns. The Thomas algorithm givesand back substitution yieldsAll three values equal $u(x_i)$ exactly: the centered finite difference for the second derivative is exact on a quadratic polynomial. This special nodal exactness does not assert exactness for arbitrary forcing or for the piecewise-linear curve joining the nodes.

For any $n$ interior points with $h=1/(n+1)$, the scaled tridiagonal matrix has $b_i=2$ and $a_i=c_i=-1$. Its pivot elements are $p_i=(i+1)/i$: start at $p_1=2$ and use $p_i=2-1/p_{i-1}$. All are positive. Interior rows have equality in diagonal dominance, so strict diagonal dominance is not the reason this example works. Instead it is a symmetric positive-definite matrix, since for a real vector $v\ne 0$ with $v_0=v_{n+1}=0$,Zero would force every neighbouring value equal and hence all values zero. This connects the boundary constraints directly to the safe pivots for tridiagonal elimination.

The displayed Thomas algorithm requires every modified pivot element to be nonzero in exact arithmetic. Nonzero original diagonal entries alone are insufficient: later pivot elements can vanish after elimination. Two useful sufficient contexts are a real symmetric positive-definite matrix, and strict diagonal dominance by rows:For the symmetric positive-definite matrix case, elimination leaves positive-definite remaining blocks by the Schur complement, so each $p_i>0$. Equivalently, the initial principal minors ${\Delta }_i=\det A_{1:i,1:i}$ are positive, and $p_i={\Delta }_i/{\Delta }_{i-1}$ with ${\Delta }_0=1$.

For strict diagonal dominance, a short induction shows $|p_i|>|c_i|$. The first row gives $|p_1|=|b_1|>|c_1|$. Assuming the preceding bound, the triangle inequality givesIn particular $p_n\ne 0$. These are structural guarantees against a zero pivot element, not universal bounds on forward error in floating-point arithmetic. A poorly conditioned linear system can still amplify roundoff error; numerical stability and the condition number must be distinguished.

An invertible matrix need not admit the unpivoted Thomas algorithm. For example,Here $\det A=-1$ and $Ax=d$, but $p_1=0$ makes $m_2=1/0$ undefined. Interchanging the two equations instead gives the upper triangular matrix system $x_1+x_2=3$, $x_2=2$, solved immediately by back substitution. Even a nonzero first diagonal does not guarantee all later pivot elements: $\left(\begin{array}{ccc}1& 1& 0\\ 1& 1& 1\\ 0& 1& 1\end{array}\right)$ has determinant $-1$, yet $p_2=1-1=0$.

For a general tridiagonal matrix, use a solver with partial pivoting when these structural guarantees are unavailable or a dangerously small pivot element occurs. Adjacent row interchanges can introduce a second upper diagonal, so a pivoted solver must track that additional fill-in; it cannot blindly continue the unpivoted recurrence relations. LAPACK's DGTSV explicitly performs partial pivoting, stores this extra upper diagonal, and includes it in back substitution. Its fixed-width storage and bounded work per row retain $O(n)$ time for one right-hand side and $O(n)$ storage.^[ref]

For the displayed unpivoted Thomas algorithm with one right-hand side, count each scalar addition, subtraction, multiplication or division as one arithmetic operation. Forming $m_i,p_i,g_i$ takes five arithmetic operations per forward step. The last unknown takes one division; each remaining back substitution step takes three arithmetic operations. ThusThis is an operation-count model, excluding assignments, indexing, pivot checks and any special simplifications of the coefficients. It is not a timing benchmark. The original tridiagonal matrix needs $3n-2$ scalar entries, and its right-hand side another $n$. Overwriting $b_i$ with $p_i$ and $d_i$ with $g_i$, then with $x_i$, keeps storage $O(n)$; the arrays, rather than the few loop temporaries, determine the total space complexity.

For repeated right-hand sides, retain the LU decomposition: computing $m_i,p_i$ once costs $3(n-1)$ arithmetic operations, and each subsequent pair of forward substitution and back substitution sweeps costs $5n-4$. Hence $r$ right-hand sides cost $3(n-1)+r(5n-4)$ arithmetic operations. This is the banded counterpart of Reusing LU factors for multiple right-hand sides.