Tridiagonal matrix algorithm

The tridiagonal matrix algorithm (TDMA), also known as the Thomas algorithm , 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,$. In matrix form, this system is written as

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. Example of such matrices commonly arise from the discretization of 1D Poisson equation (e.g., the 1D diffusion problem).

Method

The first step consists of modifying the coefficients as follows, denoting the new modified coefficients with primes:

:

:

This is the forward sweep. The solution is then obtained by back substitution:

:$x\_n\; =\; d\text{'}\_n,$

:$x\_i\; =\; d\text{'}\_i\; -\; c\text{'}\_i\; x\_\{i\; +\; 1\}\; qquad\; ;\; i\; =\; n\; -\; 1,\; n\; -\; 2,\; ldots,\; 1$

Implementation in C

The following C function will solve a general tridiagonal system. Note that the index $i$ here is zero based, in other words $i\; =\; 0,\; 1,\; dots,\; n\; -\; 1$ where $n$ is the number of unknowns.

/* Fills solution into x. Warning: will modify c and d! */void TridiagonalSolve(const double *a, const double *b, double *c, double *d, double *x, unsigned int n){ int i;

/* Modify the coefficients. */ c [0] /= b [0] ; /* Division by zero risk. */ d [0] /= b [0] ; for(i = 1; i < n; i++){ double id = 1.0/(b [i] - c [i - 1] *a [i] ); /* Division by zero risk. */ c [i] *= id; /* Last value calculated is redundant. */ d [i] = (d [i] - d [i - 1] *a [i] )*id; }

/* Now back substitute. */ x [n - 1] = d [n - 1] ; for(i = n - 2; i >= 0; i--) x [i] = d [i] - c [i] *x [i + 1] ;}

Implementation in C#

double [] TDMASolve(double [] a, double [] b, double [] c, double [] d) { double [] cc = new double [c.Length] ; double [] dd = new double [d.Length] ;

double [] x = new double [f.Length] ; int N = d.Length; cc [0] = c [0] / b [0] ; dd [0] = d [0] / b [0] ; for (int i = 1; i < N; i++) { double m = b [i] - cc [i - 1] * a [i] ; cc [i] = c [i] / m; dd [i] = (d [i] - dd [i - 1] * a [i] ) / m; } x [N - 1] = d [N - 1] ; for (int i = N - 2; i >= 0; i--) { x [i] = dd [i - 1] - cc [i - 1] * x [i + 1] ; } return x; }

Derivation

The derivation of the tridiagonal matrix algorithm involves manually performing some specialized Gaussian elimination in a generic manner.

Suppose that the unknowns are $x\_1,ldots,\; x\_n$, and that the equations to be solved are:

:

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,$

and 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:

:$ilde\; a\_i\; =\; 0,$

:$ilde\; b\_1\; =\; b\_1,$:$ilde\; b\_i\; =\; b\_i\; ilde\; b\_\{i\; -\; 1\}\; -\; ilde\; c\_\{i\; -\; 1\}\; a\_i,$

:$ilde\; c\_1\; =\; c\_1,$:$ilde\; c\_i\; =\; c\_i\; ilde\; b\_\{i\; -\; 1\},$

:$ilde\; d\_1\; =\; d\_1,$:$ilde\; d\_i\; =\; d\_i\; ilde\; b\_\{i\; -\; 1\}\; -\; ilde\; d\_\{i\; -\; 1\}\; a\_i,$

To further hasten the solution process, $ilde\; b\_i$ may be divided out (if there's no division by zero risk), the newer modified coefficients notated with an asterisk will be:

:$a\text{'}\_i\; =\; 0,$

:$b\text{'}\_i\; =\; 1,$

:$c\text{'}\_1\; =\; frac\{c\_1\}\{b\_1\},$:$c\text{'}\_i\; =\; frac\{c\_i\}\{b\_i\; -\; c\text{'}\_\{i\; -\; 1\}\; a\_i\},$

:$d\text{'}\_1\; =\; frac\{d\_1\}\{b\_1\},$:$d\text{'}\_i\; =\; frac\{d\_i\; -\; d\text{'}\_\{i\; -\; 1\}\; a\_i\}\{b\_i\; -\; c\text{'}\_\{i\; -\; 1\}\; a\_i\},$

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

:

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\text{'}\_n,$

:$x\_i\; =\; d\text{'}\_i\; -\; c\text{'}\_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:

:

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.

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.

A variant of the Thomas algorithm is the Hu-Gallash algorithm, which uses forward substitution instead of backwards substitution.

References

*cite book|author=Conte, S.D., and deBoor, C.|year=1972|title=Elementary Numerical Analysis|publisher= McGraw-Hill, New York.
*CFDWiki|name=Tridiagonal matrix algorithm &mdash; TDMA (Thomas algorithm)

External links

* [http://www.arbredeslemuriens.com/Categorie.php?IDCategorie=AlgoScilab&IDTitre=182 Thomas algorithm &mdash; SCILAB]
* [http://www.arbredeslemuriens.com/Categorie.php?IDCategorie=AlgoScilab&IDTitre=185 Thomas algorithm for tridiagonal superior periodic matrix &mdash; SCILAB]