Tridiagonal Solvers: Thomas Algorithm
Derivation and implementation of the Thomas Algorithm
The Thomas Algorithm (Thomas 1949), also known as the Tridiagonal Matrix Algorithm (TDMA), is a simplified form of Gaussian elimination specifically designed to solve tridiagonal linear systems of equations. Its efficiency stems from exploiting the sparse structure of tridiagonal matrices, leading to a much faster solution than general Gaussian elimination, \(O(N)\) operations compared to \(O(N^3)\). It is a direct method, meaning it computes the solution in a finite number of steps. The algorithm consists of two main phases: a forward elimination (or forward sweep) and a backward substitution (or backward sweep).
This article is intended for readers who are looking for a deep dive into TDMA. For an introduction to the problem, benchmarks, and discussion of algorithm choices, please see the Tridiagonal Solvers article.
Derivation of the Algorithm
The Thomas algorithm is a direct application of Gaussian elimination that transforms the tridiagonal matrix \(A\) into an upper bidiagonal matrix. The upper bidiagonal system is then easily solved by backward substitution. This process is inherently sequential and thus less suitable for parallelization than methods like Cyclic Reduction.
For this derivation, consider again our system.
\[ a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i \tag{1}\]
with the boundary coefficients \(a_0 = 0\) and \(c_{N-1} = 0\).
where: - \(b_i\) are the elements of the main diagonal. - \(a_i\) are the elements of the sub-diagonal (\(a_0\) is usually defined as 0 or ignored). - \(c_i\) are the elements of the super-diagonal (\(c_{N-1}\) is usually defined as 0 or ignored).
Forward Elimination Phase
The objective of this phase is to eliminate the sub-diagonal elements \(a_i\). We can do this by expressing each \(x_i\) as a function of only its next neighbor \(x_{i+1}\). Specifically, we aim to transform each row into a normalized form where the coefficient of \(x_i\) is 1:
\[ x_i + c'_i x_{i+1} = d'_i \tag{2}\]
For the first row \(i=0\), we know \(a_0 = 0\) and can plug in those values into Equation 1 to give
\[ b_0 x_0 + c_0 x_1 = d_0 \tag{3}\]
Dividing by \(b_0\) gives the desired form:
\[ x_0 + \left( \frac{c_0}{b_0} \right) x_1 = \frac{d_0}{b_0} \tag{4}\]
Comparing this to our goal, we define the initial modified coefficients:
\[ \begin{aligned} c'_0 &= \frac{c_0}{b_0} \\ d'_0 &= \frac{d_0}{b_0} \end{aligned} \tag{5}\]
For the remaining rows where \(i > 0\) we define a recursive process using the results from the row \(i-1\). To eliminate \(x_{i-1}\) from the \(i\)-th row, we use Equation 2 for the previous row \((i-1)\):
\[ x_{i-1} + c'_{i-1} x_i = d'_{i-1} \tag{6}\]
Rearranging this gives an expression for \(x_{i-1}\) in terms of \(x_i\):
\[ x_{i-1} = d'_{i-1} - c'_{i-1} x_i \tag{7}\]
We can now substitute this expression into the general tridiagonal equation Equation 1 for row \(i\):
\[ a_i (d'_{i-1} - c'_{i-1} x_i) + b_i x_i + c_i x_{i+1} = d_i \tag{8}\]
Rearrange to group terms with \(x_i\):
\[ (b_i - a_i c'_{i-1}) x_i + c_i x_{i+1} = d_i - a_i d'_{i-1} \tag{9}\]
To reach our goal form (\(x_i + c'_i x_{i+1} = d'_i\)), we divide by the coefficient \((b_i - a_i c'_{i-1})\):
\[ x_i + \left( \frac{c_i}{b_i - a_i c'_{i-1}} \right) x_{i+1} = \frac{d_i - a_i d'_{i-1}}{b_i - a_i c'_{i-1}} \tag{10}\]
Thus, for \(i = 1, \dots, N-1\), the recurrence relations for the modified coefficients are:
\[ \begin{aligned} c'_i &= \frac{c_i}{b_i - a_i c'_{i-1}} \\ d'_i &= \frac{d_i - a_i d'_{i-1}}{b_i - a_i c'_{i-1}} \end{aligned} \tag{11}\]
The term \((b_i - a_i c'_{i-1})\) is known as the pivot. It must be non-zero for the algorithm to succeed. For most physical problems, strict diagonal dominance ensures stability.
Backward Substitution Phase
Once the forward elimination is complete, the system is in an upper bidiagonal form:
\[ x_i + c'_i x_{i+1} = d'_i \tag{12}\]
Starting at the last row (\(i = N-1\)), since \(c_{N-1} = 0\), we have \(c'_{N-1} = 0\). The solution for the last unknown is simply:
\[ x_{N-1} = d'_{N-1} \tag{13}\]
We can then solve for all remaining unknowns by working backward from \(i = N-2\) down to \(0\):
\[ x_i = d'_i - c'_i x_{i+1} \tag{14}\]
By computing \(x_{N-1}\) first and then iteratively finding \(x_{N-2}, x_{N-3}, ..., x_0\), we obtain the complete solution.
Pseudocode for Implementation
function solve_tridiag_thomas_algorithm(a, b, c, d):
n = length(d)
// Create temporary arrays for modified coefficients
cp = array_of_zeros(n - 1) // c'
dp = array_of_zeros(n) // d'
x = array_of_zeros(n) // Solution vector
// --- Phase 1: Forward Elimination ---
// Handle the first equation (i=0)
cp[0] = c[0] / b[0]
dp[0] = d[0] / b[0]
// Iterate for i = 1 to N-1
for i from 1 to n - 1:
denom = b[i] - a[i-1] * cp[i-1]
// If denom is zero, the system is singular or ill-conditioned
if i < n - 1: // Only calculate cp[i] if c[i] exists
cp[i] = c[i] / denom
dp[i] = (d[i] - a[i-1] * dp[i-1]) / denom
// --- Phase 2: Backward Substitution ---
// Solve for the last unknown (x_{N-1})
x[n-1] = dp[n-1]
// Iterate for i = N-2 down to 0
for i from n - 2 down to 0:
x[i] = dp[i] - cp[i] * x[i+1]
return x