Thomas algorithm

The Thomas algorithm (after Llewellyn Thomas ) or also tridiagonal matrix algorithm (TDMA) is a simplified form of the Gaussian elimination method , which is used to quickly solve linear systems of equations with a tridiagonal matrix .

Procedure

A tridiagonal system with n unknowns can be written as ${\displaystyle x_{i}\,\!}$

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

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

In matrix form the system is written as: ${\displaystyle \left[{\begin{matrix}{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{matrix}}\right]\left[{\begin{matrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\\vdots \\{x_{n}}\\\end{matrix}}\right]=\left[{\begin{matrix}{d_{1}}\\{d_{2}}\\{d_{3}}\\\vdots \\{d_{n}}\\\end{matrix}}\right].}$

Examples of these matrices usually arise from the discretization of the one-dimensional Poisson equation (e.g. one-dimensional diffusion problems ) and from natural cubic spline interpolation .

For such systems, the solution can be found with the Thomas algorithm after operations: A first pass eliminates the s, then the solution is obtained with the help of a backward insertion process . ${\displaystyle {\mathcal {O}}(n)}$ ${\displaystyle a_{i}}$

The first step is to recursively modify the coefficients as follows (the "new" modified coefficients are marked with a dash):

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

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

This is the forward pass. The solution then results from a backward insertion process:

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

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

variants

In some situations, especially in those with periodic boundary conditions , it may be that a slightly different form of the tridiagonal system has to be solved:

${\displaystyle {\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 matrix form: ${\displaystyle \left[{\begin{matrix}{b_{1}}&{c_{1}}&{}&{}&{a_{1}}\\{a_{2}}&{b_{2}}&{c_{2}}&{}&{}\\{}&{a_{3}}&{b_{3}}&\ddots &{}\\{}&{}&\ddots &\ddots &{c_{n-1}}\\{c_{n}}&{}&{}&{a_{n}}&{b_{n}}\\\end{matrix}}\right]\left[{\begin{matrix}{x_{1}}\\{x_{2}}\\{x_{3}}\\\vdots \\{x_{n}}\\\end{matrix}}\right]=\left[{\begin{matrix}{d_{1}}\\{d_{2}}\\{d_{3}}\\\vdots \\{d_{n}}\\\end{matrix}}\right].}$

In this case the Sherman-Morrison-Woodbury formula can be used to utilize the Thomas algorithm and avoid the additional operations of the Gaussian elimination method.

In other cases the system of equations can be block- tridiagonal, i.e. H. the elements of the above system of equations are smaller block matrices (e.g. in the two-dimensional Poisson equation ). Simplified variants of Gaussian elimination have also been developed for these cases.

Another variant of the Thomas algorithm is the Hu-Gallash algorithm, which uses a forward insertion method instead of the backward insertion method.

Derivation

Let the unknowns be . The equations to be solved are: ${\displaystyle x_{1},\ldots ,x_{n}}$

${\displaystyle {\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}}}$

The second equation ( ) is now modified with the first equation as follows: ${\displaystyle i=2}$

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

This gives:

${\displaystyle {\begin{aligned}(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}\,\end{aligned}}}$

This eliminated from the second equation. Using an analogous procedure, the modified second and third equations result: ${\displaystyle x_{1}}$

${\displaystyle {\begin{aligned}(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}\,\end{aligned}}}$

Now has been eliminated. If this procedure is repeated up to the -th line, the modified -th equation contains only one unknown. This equation is now solved and the result used to solve the -th equation. This continues until all unknowns are known. ${\displaystyle x_{2}}$ ${\displaystyle n}$ ${\displaystyle n}$ ${\displaystyle (n-1)}$

Understandably, the coefficients get more complicated with each step when they are explicitly represented. However, the coefficients can also be represented recursively:

${\displaystyle {\begin{aligned}{\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}\,\end{aligned}}}$

In order to further accelerate the solution process, it is divided out (if ). The new coefficients are: ${\displaystyle {\tilde {b}}_{i}}$ ${\displaystyle {\tilde {b}}_{i}\neq 0}$

${\displaystyle {\begin{aligned}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}}}\,\end{aligned}}}$

This gives:

${\displaystyle {\begin{array}{llcl}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 contains only one unknown. Finding them reduces the number of unknowns in the penultimate equation to one so that this backward insertion method can be used to determine all of the unknowns.

${\displaystyle {\begin{aligned}x_{n}&=d'_{n}\\x_{i}&=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1\end{aligned}}}$

Individual evidence

• Conte, SD, and deBoor, C .: Elementary Numerical Analysis . McGraw-Hill, New York., 1972.