Numerical linear algebra is the study of algorithms for performing linear algebra computations, most notably matrix operations, on computers. It is often a fundamental part of engineering and computational science problems, such as image and signal processing, Telecommunication, computational...
, the tridiagonal matrix algorithm (TDMA), also known as the Thomas algorithm (named after Llewellyn Thomas
Llewellyn Thomas
Llewellyn Hilleth Thomas was a British physicist and applied mathematician. He is best known for his contributions to atomic physics, in particular:...
In linear algebra, Gaussian elimination is an algorithm for solving systems of linear equations. It can also be used to find the rank of a matrix, to calculate the determinant of a matrix, and to calculate the inverse of an invertible square matrix...
that can be used to solve tridiagonal systems of equations. A tridiagonal system for n unknowns may be written as
where and . In matrix form, this system is written as
For such systems, the solution can be obtained in operations instead of required by Gaussian elimination
Gaussian elimination
In linear algebra, Gaussian elimination is an algorithm for solving systems of linear equations. It can also be used to find the rank of a matrix, to calculate the determinant of a matrix, and to calculate the inverse of an invertible square matrix...
. A first sweep eliminates the '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
Heat equation
The heat equation is an important partial differential equation which describes the distribution of heat in a given region over time...
In the mathematical field of numerical analysis, spline interpolation is a form of interpolation where the interpolant is a special type of piecewise polynomial called a spline. Spline interpolation is preferred over polynomial interpolation because the interpolation error can be made small even...
.
Method
The first step consists of modifying the coefficients as follows, denoting the new modified coefficients with primes:
and:
This is the forward sweep. The solution is then obtained by back substitution:
C is a general-purpose computer programming language developed between 1969 and 1973 by Dennis Ritchie at the Bell Telephone Laboratories for use with the Unix operating system....
function will solve a general tridiagonal system (though it will destroy the input vectors b and v in the process). Note that the index here is zero based
Array data type
In computer science, an array type is a data type that is meant to describe a collection of elements , each selected by one or more indices that can be computed at run time by the program. Such a collection is usually called an array variable, array value, or simply array...
, in other words where is the number of unknowns.
void solveMatrix (int n, double *a, double *b, double *c, double *v, double *x)
{
/**
* n - number of equations
* a - sub-diagonal (means it is the diagonal below the main diagonal) -- indexed from 1..n-1
* b - the main diagonal
* c - sup-diagonal (means it is the diagonal above the main diagonal) -- indexed from 0..n-2
* v - right part
* x - the answer
*/
for (int i = 1; i < n; i++)
{
double m = a[i]/b[i-1];
b[i] = b[i] - m*c[i-1];
v[i] = v[i] - m*v[i-1];
}
x[n-1] = v[n-1]/b[n-1];
for (int i = n - 2; i >= 0; i--)
x[i]=(v[i]-c[i]*x[i+1])/b[i];
}
In computer science, an array type is a data type that is meant to describe a collection of elements , each selected by one or more indices that can be computed at run time by the program. Such a collection is usually called an array variable, array value, or simply array...
, in other words where is the number of unknowns.
function x = TDMAsolver(a,b,c,d)
%a, b, c are the column vectors for the compressed tridiagonal matrix, d is the right vector
n = length(b); % n is the number of rows
% Modify the first-row coefficients
c(1) = c(1) / b(1); % Division by zero risk.
d(1) = d(1) / b(1); % Division by zero would imply a singular matrix.
for i = 2:n
id = 1 / (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) - d(i-1) * a(i)) * id;
end
% Now back substitute.
x(n) = d(n);
for i = n-1:-1:1
x(i) = d(i) - c(i) * x(i + 1);
end
In computer science, an array type is a data type that is meant to describe a collection of elements , each selected by one or more indices that can be computed at run time by the program. Such a collection is usually called an array variable, array value, or simply array...
, in other words where is the number of unknowns.
Sometimes it is undesirable to have the solver routine overwrite the tridiagonal coefficients (e.g. for solving multiple systems of equations where only the right side of the system changes), so this implementation gives an example of a relatively inexpensive method of preserving the coefficients.
subroutine solve_tridiag(a,b,c,v,x,n)
implicit none
! a - sub-diagonal (means it is the diagonal below the main diagonal)
! b - the main diagonal
! c - sup-diagonal (means it is the diagonal above the main diagonal)
! v - right part
! x - the answer
! n - number of equations
integer,intent(in) :: n
real(8),dimension(n),intent(in) :: a,b,c,v
real(8),dimension(n),intent(out) :: x
real(8),dimension(n) :: bp,vp
real(8) :: m
integer i
! Make copies of the b and v variables so that they are unaltered by this sub
bp(1) = b(1)
vp(1) = v(1)
!The first pass (setting coefficients):
firstpass: do i = 2,n
m = a(i)/bp(i-1)
bp(i) = b(i) - m*c(i-1)
vp(i) = v(i) - m*vp(i-1)
end do firstpass
x(n) = vp(n)/bp(n)
!The second pass (back-substition)
backsub:do i = n-1, 1, -1
x(i) = (vp(i) - c(i)*x(i+1))/bp(i)
end do backsub
end subroutine solve_tridiag
Derivation
The derivation of the tridiagonal matrix algorithm involves manually performing some specialized Gaussian elimination
Gaussian elimination
In linear algebra, Gaussian elimination is an algorithm for solving systems of linear equations. It can also be used to find the rank of a matrix, to calculate the determinant of a matrix, and to calculate the inverse of an invertible square matrix...
in a generic manner.
Suppose that the unknowns are , and that the equations to be solved are:
Consider modifying the second () equation with the first equation as follows:
which would give:
and the effect is that has been eliminated from the second equation. Using a similar tactic with the modified second equation on the third equation yields:
This time was eliminated. If this procedure is repeated until the row; the (modified) equation will involve only one unknown, . This may be solved for and then used to solve the 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:
To further hasten the solution process, may be divided out (if there's no division by zero risk), the newer modified coefficients notated with an asterisk will be:
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:
Variants
In some situations, particularly those involving periodic boundary conditions, a slightly perturbed form of the tridiagonal system may need to be solved:
In mathematics, in particular linear algebra, the Sherman–Morrison formula, named after Jack Sherman and Winifred J. Morrison, computes the inverse of the sum of an invertiblematrix Aand the dyadic product, u v^T,of a column vector u and a row vector v^T...
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
Block matrix
In the mathematical discipline of matrix theory, a block matrix or a partitioned matrix is a matrix broken into sections called blocks. Looking at it another way, the matrix is written in terms of smaller matrices. We group the rows and columns into adjacent 'bunches'. A partition is the rectangle...
), 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.
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.
and 
. In matrix form, this system is written as
operations instead of 
required by Gaussian elimination
'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
here is zero based
where 
is the number of unknowns.
here is one based
where 
is the number of unknowns.
here is one based
where 
is the number of unknowns.
, and that the equations to be solved are:
) equation with the first equation as follows:
has been eliminated from the second equation. Using a similar tactic with the modified second equation on the third equation yields:
was eliminated. If this procedure is repeated until the 
row; the (modified) 
equation will involve only one unknown, 
. This may be solved for and then used to solve the 
equation, and so on until all of the unknowns are solved for.
may be divided out (if there's no division by zero risk), the newer modified coefficients notated with an asterisk will be:
