Tridiagonal matrix algorithm

From Wikipedia, the free encyclopedia

Jump to navigationJump to search

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

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

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

${\displaystyle {\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 ${\displaystyle O(n)}$ operations instead of ${\displaystyle O(n^{3})}$ required by Gaussian elimination. A first sweep eliminates the ${\displaystyle 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 and natural cubic spline interpolation.

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]

Contents

• 1Method
• 2Derivation
• 3Variants
• 4References

Method[edit]

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

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

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

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

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

The method above preserves the original coefficient vectors. If this is not required, then a much simpler form of the algorithm is

${\displaystyle {\begin{array}{lcl}w_{i}={\cfrac {a_{i}}{b_{i-1}}}\\b_{i}:=b_{i}-w_{i}c_{i-1}\\d_{i}:=d_{i}-w_{i}d_{i-1}\end{array}},i=2,3,\dots ,n}$

followed by the back substitution

${\displaystyle {\begin{array}{lcl}x_{n}={\cfrac {d_{n}}{b_{n}}}\\x_{i}={\cfrac {d_{i}-c_{i}x_{i+1}}{b_{i}}},i=n-1,n-2,\dots ,1\end{array}}}$

The implementation in a VBA subroutine without preserving the coefficient vectors is shown below

Sub TriDiagonal_Matrix_Algorithm(N%, A#(), B#(), C#(), D#(), X#())
 Dim i%, W#
 For i = 2 To N
  W = A(i) / B(i - 1)
  B(i) = B(i) - W * C(i - 1)
  D(i) = D(i) - W * D(i - 1)
 Next i
 X(N) = D(N) / B(N)
 For i = N - 1 To 1 Step -1
  X(i) = (D(i) - C(i) * X(i + 1)) / B(i)
 Next i
End Sub