Tridiagonal matrix algorithm

In numerical linear algebra, the tridiagonal matrix algorithm, also known as the Thomas algorithm (named after Llewellyn Thomas), is a simplified form of Gaussian elimination that can be used to solve tridiagonal systems of equations. A tridiagonal system for n unknowns may be written as

$a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},\,\!$

where $a_{1}=0\,$ and $c_{n}=0\,$.

${\begin{bmatrix}{b_{1}}&{c_{1}}&{}&{}&{0}\\{a_{2}}&{b_{2}}&{c_{2}}&{}&{}\\{}&{a_{3}}&{b_{3}}&\ddots &{}\\{}&{}&\ddots &\ddots &{c_{n-1}}\\{0}&{}&{}&{a_{n}}&{b_{n}}\\\end{bmatrix}}{\begin{bmatrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\\vdots \\{x_{n}}\\\end{bmatrix}}={\begin{bmatrix}{d_{1}}\\{d_{2}}\\{d_{3}}\\\vdots \\{d_{n}}\\\end{bmatrix}}.$

For such systems, the solution can be obtained in $O(n)$ operations instead of $O(n^{3})$ required by Gaussian elimination. A first sweep eliminates the $a_{i}$'s, and then an (abbreviated) backward substitution produces the solution. Examples of such matrices commonly arise from the discretization of 1D Poisson equation (e.g., the 1D diffusion problem) and natural cubic spline interpolation; similar systems of matrices arise in tight binding physics or nearest neighbor effects models.

Thomas' algorithm is not stable in general, but is so in several special cases, such as when the matrix is diagonally dominant (either by rows or columns) or symmetric positive definite;^[1][2] for a more precise characterization of stability of Thomas' algorithm, see Higham Theorem 9.12.^[3] If stability is required in the general case, Gaussian elimination with partial pivoting (GEPP) is recommended instead.^[2]

Method

The forward sweep consists of modifying the coefficients as follows, denoting the new coefficients with primes:

$c'_{i}={\begin{cases}{\begin{array}{lcl}{\cfrac {c_{i}}{b_{i}}}&;&i=1\\{\cfrac {c_{i}}{b_{i}-a_{i}c'_{i-1}}}&;&i=2,3,\dots ,n-1\\\end{array}}\end{cases}}\,$

and

$d'_{i}={\begin{cases}{\begin{array}{lcl}{\cfrac {d_{i}}{b_{i}}}&;&i=1\\{\cfrac {d_{i}-a_{i}d'_{i-1}}{b_{i}-a_{i}c'_{i-1}}}&;&i=2,3,\dots ,n.\\\end{array}}\end{cases}}\,$

The solution is then obtained by back substitution:

$x_{n}=d'_{n}\,$

$x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1.$

Derivation

The derivation of the tridiagonal matrix algorithm is a special case of Gaussian elimination.

Suppose that the unknowns are $x_{1},\ldots ,x_{n}$, and that the equations to be solved are:

${\begin{aligned}b_{1}x_{1}+c_{1}x_{2}&=d_{1};&i&=1\\a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}&=d_{i};&i&=2,\ldots ,n-1\\a_{n}x_{n-1}+b_{n}x_{n}&=d_{n};&i&=n.\end{aligned}}$

Consider modifying the second ($i=2$) equation with the first equation as follows:

$({\mbox{equation 2}})\cdot b_{1}-({\mbox{equation 1}})\cdot a_{2}$

which would give:

$(a_{2}x_{1}+b_{2}x_{2}+c_{2}x_{3})b_{1}-(b_{1}x_{1}+c_{1}x_{2})a_{2}=d_{2}b_{1}-d_{1}a_{2}\,$

$(b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3}=d_{2}b_{1}-d_{1}a_{2}\,$

where the second equation immediately above is a simplified version of the equation immediately preceding it. The effect is that $x_{1}$ has been eliminated from the second equation. Using a similar tactic with the modified second equation on the third equation yields:

$(a_{3}x_{2}+b_{3}x_{3}+c_{3}x_{4})(b_{2}b_{1}-c_{1}a_{2})-((b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3})a_{3}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}\,$

$(b_{3}(b_{2}b_{1}-c_{1}a_{2})-c_{2}b_{1}a_{3})x_{3}+c_{3}(b_{2}b_{1}-c_{1}a_{2})x_{4}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}.\,$

This time $x_{2}$ was eliminated. If this procedure is repeated until the $n^{th}$ row; the (modified) $n^{th}$ equation will involve only one unknown, $x_{n}$. This may be solved for and then used to solve the $(n-1)^{th}$ equation, and so on until all of the unknowns are solved for.

Clearly, the coefficients on the modified equations get more and more complicated if stated explicitly. By examining the procedure, the modified coefficients (notated with tildes) may instead be defined recursively:

${\tilde {a}}_{i}=0\,$

${\tilde {b}}_{1}=b_{1}\,$

${\tilde {b}}_{i}=b_{i}{\tilde {b}}_{i-1}-{\tilde {c}}_{i-1}a_{i}\,$

${\tilde {c}}_{1}=c_{1}\,$

${\tilde {c}}_{i}=c_{i}{\tilde {b}}_{i-1}\,$

${\tilde {d}}_{1}=d_{1}\,$

${\tilde {d}}_{i}=d_{i}{\tilde {b}}_{i-1}-{\tilde {d}}_{i-1}a_{i}.\,$

To further hasten the solution process, ${\tilde {b}}_{i}$ may be divided out (if there's no division by zero risk), the newer modified coefficients, each notated with a prime, will be:

$a'_{i}=0\,$

$b'_{i}=1\,$

$c'_{1}={\frac {c_{1}}{b_{1}}}\,$

$c'_{i}={\frac {c_{i}}{b_{i}-c'_{i-1}a_{i}}}\,$

$d'_{1}={\frac {d_{1}}{b_{1}}}\,$

$d'_{i}={\frac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}}.\,$

This gives the following system with the same unknowns and coefficients defined in terms of the original ones above:

${\begin{array}{lcl}x_{i}+c'_{i}x_{i+1}=d'_{i}\qquad &;&\ i=1,\ldots ,n-1\\x_{n}=d'_{n}\qquad &;&\ i=n.\\\end{array}}\,$

The last equation involves only one unknown. Solving it in turn reduces the next last equation to one unknown, so that this backward substitution can be used to find all of the unknowns:

$x_{n}=d'_{n}\,$

$x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1.$

Variants

In some situations, particularly those involving periodic boundary conditions, a slightly perturbed form of the tridiagonal system may need to be solved:

${\begin{aligned}a_{1}x_{n}+b_{1}x_{1}+c_{1}x_{2}&=d_{1},\\a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}&=d_{i},\quad \quad i=2,\ldots ,n-1\\a_{n}x_{n-1}+b_{n}x_{n}+c_{n}x_{1}&=d_{n}.\end{aligned}}$

In this case, we can make use of the Sherman-Morrison formula to avoid the additional operations of Gaussian elimination and still use the Thomas algorithm. The method requires solving a modified non-cyclic version of the system for both the input and a sparse corrective vector, and then combining the solutions. This can be done efficiently if both solutions are computed at once, as the forward portion of the pure tridiagonal matrix algorithm can be shared.

In other situations, the system of equations may be block tridiagonal (see block matrix), with smaller submatrices arranged as the individual elements in the above matrix system(e.g., the 2D Poisson problem). Simplified forms of Gaussian elimination have been developed for these situations.^[4]

The textbook Numerical Mathematics by Quarteroni, Sacco and Saleri, lists a modified version of the algorithm which avoids some of the divisions (using instead multiplications), which is beneficial on some computer architectures.

References

1. Pradip Niyogi (2006). Introduction to Computational Fluid Dynamics. Pearson Education India. p. 76. ISBN 978-81-7758-764-7.
2. Biswa Nath Datta (2010). Numerical Linear Algebra and Applications, Second Edition. SIAM. p. 162. ISBN 978-0-89871-765-5.
3. Nicholas J. Higham (2002). Accuracy and Stability of Numerical Algorithms: Second Edition. SIAM. p. 175. ISBN 978-0-89871-802-7.
4. Quarteroni, Alfio; Sacco, Riccardo; Saleri, Fausto (2007). "Section 3.8". Numerical Mathematics. Springer, New York. ISBN 978-3-540-34658-6.

• Conte, S.D. & deBoor, C. (1972). Elementary Numerical Analysis. McGraw-Hill, New York. ISBN 0070124469.
• This article incorporates text from the article Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) on CFD-Wiki that is under the GFDL license.

• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). "Section 2.4". Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8