Tridiagonal matrix algorithm

From Wikipedia, the free encyclopedia

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

$\left[ \begin{matrix} {b_1} & {c_1} & { } & { } & { 0 } \\ {a_2} & {b_2} & {c_2} & { } & { } \\ { } & {a_3} & {b_3} & \cdot & { } \\ { } & { } & \cdot & \cdot & {c_{n-1}}\\ { 0 } & { } & { } & {a_n} & {b_n}\\ \end{matrix} \right] \left[ \begin{matrix} {x_1 } \\ {x_2 } \\ \cdot \\ \cdot \\ {x_n } \\ \end{matrix} \right] = \left[ \begin{matrix} {d_1 } \\ {d_2 } \\ \cdot \\ \cdot \\ {d_n } \\ \end{matrix} \right].$

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).

1 Method
2 Variants
3 References
4 External links

[edit] Method

See the derivation.

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

$a'_i = 0\,$

$b'_i = 1\,$

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

$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 sweep. The solution is then obtained by back substitution:

$x_n = d'_n\,$

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

[edit] 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] = c[0]/b[0];                              /* Division by zero risk. */
        d[0] = d[0]/b[0];
        double id;
        for(i = 1; i != n; i++){
                id = 1.0/(b[i] - c[i - 1]*a[i]);      /* Division by zero risk. */
                c[i] = c[i]*id;                             /* Last value calculated is redundant. */
                d[i] = (d[i] - a[i]*d[i - 1])*id;
        }
 
        /* Now back substitute. */
        x[n - 1] = d[n - 1];
        for(i = n - 2; i != -1; i--)
                x[i] = d[i] - c[i]*x[i + 1];
}

[edit] 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];
            double m = 0D;
            for (int i = 1; i < N; i++)
            {
                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;
        }

[edit] Implementation in Matlab

function [x]=TDMAsolver(a,b,c,d)
    n=length(a);
 
    c(1) = c(1)/b(1);
    d(1) = d(1)/b(1);
    for i = 2:n
        id = (b(i) - c(i-1)*a(i));
        c(i)=c(i)/id; 
        d(i)=(d(i) - a(i)*d(i-1))/id;
    end
 
    x(n) = d(n);
    for i= n-1:-1:1
        x(i) = d(i) - c(i)*x(i + 1);
    end

[edit] Variants

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

$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. \,\!$

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.

[edit] References

Conte, S.D., and deBoor, C. (1972). Elementary Numerical Analysis. McGraw-Hill, New York..
This article incorporates text from the article matrix algorithm - TDMA (Thomas algorithm) Tridiagonal matrix algorithm - TDMA (Thomas algorithm) on CFD-Wiki that is under the GFDL license.