|
|
| Line 1: |
Line 1: |
| 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 matrix|tridiagonal systems of equations]]. A tridiagonal system for ''n'' unknowns may be written as
| | Hello. Allow me introduce the author. Her title is Refugia Shryock. To perform baseball is the pastime he will by no means stop performing. Managing individuals is his profession. California is our birth location.<br><br>Have a look at my blog [http://www.noonptm.com/fox/upload/index.php?do=/profile-3715/info/ over the counter std test] |
| | |
| :<math>
| |
| a_i x_{i - 1} + b_i x_i + c_i x_{i + 1} = d_i , \,\!</math>
| |
| where <math> a_1 = 0\, </math> and <math> c_n = 0\, </math>.
| |
| | |
| :<math>
| |
| \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}
| |
| .
| |
| </math>
| |
| | |
| For such systems, the solution can be obtained in <math>O(n)</math> operations instead of <math>O(n^3)</math> required by [[Gaussian elimination]]. A first sweep eliminates the <math>a_i</math>'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 [[Heat equation|diffusion problem]]) and natural cubic [[spline interpolation]]; similar systems of matrices arise in [[tight binding |tight binding physics]] or [[nearest neighbor]] effects models.
| |
| | |
| ==Method==
| |
| The forward sweep consists of modifying the coefficients as follows, denoting the new modified coefficients with primes:
| |
| | |
| :<math>c'_i =
| |
| \begin{cases}
| |
| \begin{array}{lcl}
| |
| \cfrac{c_i}{b_i} & ; & i = 1 \\
| |
| \cfrac{c_i}{b_i - c'_{i - 1} a_i} & ; & i = 2, 3, \dots, n-1 \\
| |
| \end{array}
| |
| \end{cases}
| |
| \,</math>
| |
| | |
| and
| |
| | |
| :<math>d'_i =
| |
| \begin{cases}
| |
| \begin{array}{lcl}
| |
| \cfrac{d_i}{b_i} & ; & 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}
| |
| \,</math>
| |
| | |
| The solution is then obtained by back substitution:
| |
| | |
| :<math>x_n = d'_n\,</math>
| |
| | |
| :<math>x_i = d'_i - c'_i x_{i + 1} \qquad ; \ i = n - 1, n - 2, \ldots, 1.</math>
| |
| | |
| ==Implementations==
| |
| All the provided implementations assume that the three diagonals, a (below), b (main), and c (above), are passed as arguments.
| |
| | |
| ===C===
| |
| The following [[C (programming language)|C]] function will solve a general tridiagonal system (though it will destroy the input vector c in the process). Note that the index <math>i</math> here is [[Array data type|zero based]], in other words <math>i = 0, 1, \dots, N - 1</math> where <math>N</math> is the number of unknowns.
| |
| | |
| <source lang="c" enclose="div">
| |
| void solve_tridiagonal_in_place_destructive(double x[], const size_t N, const double a[], const double b[], double c[]) {
| |
| /* unsigned integer of same size as pointer */
| |
| size_t in;
| |
|
| |
| /*
| |
| solves Ax = v where A is a tridiagonal matrix consisting of vectors a, b, c
| |
| note that contents of input vector c will be modified, making this a one-time-use function
| |
| x[] - initially contains the input vector v, and returns the solution x. indexed from [0, ..., N - 1]
| |
| N — number of equations
| |
| a[] - subdiagonal (means it is the diagonal below the main diagonal) -- indexed from [1, ..., N - 1]
| |
| b[] - the main diagonal, indexed from [0, ..., N - 1]
| |
| c[] - superdiagonal (means it is the diagonal above the main diagonal) -- indexed from [0, ..., N - 2]
| |
| */
| |
|
| |
| c[0] = c[0] / b[0];
| |
| x[0] = x[0] / b[0];
| |
|
| |
| /* loop from 1 to N - 1 inclusive */
| |
| for (in = 1; in < N; in++) {
| |
| double m = 1.0 / (b[in] - a[in] * c[in - 1]);
| |
| c[in] = c[in] * m;
| |
| x[in] = (x[in] - a[in] * x[in - 1]) * m;
| |
| }
| |
|
| |
| /* loop from N - 2 to 0 inclusive, safely testing loop end condition */
| |
| for (in = N - 1; in-- > 0; )
| |
| x[in] = x[in] - c[in] * x[in + 1];
| |
| }
| |
| | |
| </source>
| |
| | |
| The following variant preserves the system of equations for reuse on other inputs. Note the necessity of library calls to allocate and free scratch space - a more efficient implementation for solving the same tridiagonal system on many inputs would rely on the calling function to provide a pointer to the scratch space.
| |
| <source lang="c">
| |
| void solve_tridiagonal_in_place_reusable(double x[], const size_t N, const double a[], const double b[], const double c[]) {
| |
| size_t in;
| |
|
| |
| /* Allocate scratch space. */
| |
| double * cprime = malloc(sizeof(double) * N);
| |
|
| |
| if (!cprime) {
| |
| /* do something to handle error */
| |
| }
| |
|
| |
| cprime[0] = c[0] / b[0];
| |
| x[0] = x[0] / b[0];
| |
|
| |
| /* loop from 1 to N - 1 inclusive */
| |
| for (in = 1; in < N; in++) {
| |
| double m = 1.0 / (b[in] - a[in] * cprime[in - 1]);
| |
| cprime[in] = c[in] * m;
| |
| x[in] = (x[in] - a[in] * x[in - 1]) * m;
| |
| }
| |
|
| |
| /* loop from N - 2 to 0 inclusive, safely testing loop end condition */
| |
| for (in = N - 1; in-- > 0; )
| |
| x[in] = x[in] - cprime[in] * x[in + 1];
| |
|
| |
| /* free scratch space */
| |
| free(cprime);
| |
| }
| |
| </source>
| |
| | |
| ===Python===
| |
| Note that the index <math>i</math> here is [[Array data type|zero-based]], in other words <math>i = 0, 1, \dots, N-1</math> where <math>N</math> is the number of unknowns.
| |
| | |
| <source lang="Python">
| |
| # note: function also modifies b[] and d[] params while solving
| |
| def TDMASolve(a, b, c, d):
| |
| n = len(d) # n is the numbers of rows, a and c has length n-1
| |
| for i in xrange(n-1):
| |
| d[i+1] -= d[i] * a[i] / b[i]
| |
| b[i+1] -= c[i] * a[i] / b[i]
| |
| for i in reversed(xrange(n-1)):
| |
| d[i] -= d[i+1] * c[i] / b[i+1]
| |
| return [d[i] / b[i] for i in xrange(n)] # return the solution
| |
| </source>
| |
| | |
| ===MATLAB===
| |
| <source lang="MATLAB">
| |
| 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(d); % 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);
| |
|
| |
| for i = 2:n-1
| |
| temp = b(i) - a(i) * c(i-1);
| |
| c(i) = c(i) / temp;
| |
| d(i) = (d(i) - a(i) * d(i-1))/temp;
| |
| end
| |
|
| |
| d(n) = (d(n) - a(n) * d(n-1))/( b(n) - a(n) * c(n-1));
| |
|
| |
| % Now back substitute.
| |
| x(n) = d(n);
| |
| for i = n-1:-1:1
| |
| x(i) = d(i) - c(i) * x(i + 1);
| |
| end
| |
| </source>
| |
| | |
| ===Fortran 90===
| |
| Note that the index <math>i</math> here is [[Array data type|one based]], in other words <math>i = 1, 2, \dots, n</math> where <math>n</math> 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.
| |
| | |
| <source lang="Fortran">
| |
| subroutine solve_tridiag(a,b,c,d,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)
| |
| ! d - right part
| |
| ! x - the answer
| |
| ! n - number of equations
| |
| | |
| integer,intent(in) :: n
| |
| real(8),dimension(n),intent(in) :: a,b,c,d
| |
| real(8),dimension(n),intent(out) :: x
| |
| real(8),dimension(n) :: cp,dp
| |
| real(8) :: m
| |
| integer i
| |
| | |
| ! initialize c-prime and d-prime
| |
| cp(1) = c(1)/b(1)
| |
| dp(1) = d(1)/b(1)
| |
| ! solve for vectors c-prime and d-prime
| |
| do i = 2,n
| |
| m = b(i)-cp(i-1)*a(i)
| |
| cp(i) = c(i)/m
| |
| dp(i) = (d(i)-dp(i-1)*a(i))/m
| |
| enddo
| |
| ! initialize x
| |
| x(n) = dp(n)
| |
| ! solve for x from the vectors c-prime and d-prime
| |
| do i = n-1, 1, -1
| |
| x(i) = dp(i)-cp(i)*x(i+1)
| |
| end do
| |
| | |
| end subroutine solve_tridiag
| |
| </source>
| |
| | |
| This subroutine offers an option of overwriting d or not.<ref>*[http://www.ivt.ntnu.no/ept/mtf/fagtilbud/tkt4140/diskett/ /App9/tdma.f90]</ref>
| |
| | |
| <source lang="Fortran">
| |
| subroutine tdma(n,a,b,c,d,x)
| |
| implicit none
| |
| integer, intent(in) :: n
| |
| real, intent(in) :: a(n), c(n)
| |
| real, intent(inout), dimension(n) :: b, d
| |
| real, intent(out) :: x(n)
| |
| ! --- Local variables ---
| |
| integer :: i
| |
| real :: q
| |
| ! --- Elimination ---
| |
| do i = 2,n
| |
| q = a(i)/b(i - 1)
| |
| b(i) = b(i) - c(i - 1)*q
| |
| d(i) = d(i) - d(i - 1)*q
| |
| end do
| |
| ! --- Backsubstitution ---
| |
| q = d(n)/b(n)
| |
| x(n) = q
| |
| do i = n - 1,1,-1
| |
| q = (d(i) - c(i)*q)/b(i)
| |
| x(i) = q
| |
| end do
| |
| return
| |
| end
| |
| </source>
| |
| | |
| ==Derivation==
| |
| The derivation of the tridiagonal matrix algorithm involves manually performing some specialized [[Gaussian elimination]] in a generic manner.
| |
| | |
| Suppose that the unknowns are <math>x_1,\ldots, x_n</math>, and that the equations to be solved are:
| |
| | |
| :<math>\begin{align}
| |
| 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{align}
| |
| </math>
| |
| | |
| Consider modifying the second (<math>i = 2</math>) equation with the first equation as follows:
| |
| | |
| :<math>
| |
| (\mbox{equation 2}) \cdot b_1 - (\mbox{equation 1}) \cdot a_2
| |
| </math>
| |
| | |
| which would give:
| |
| | |
| :<math>
| |
| (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
| |
| \,</math>
| |
| | |
| :<math>
| |
| (b_2 b_1 - c_1 a_2) x_2 + c_2 b_1 x_3 = d_2 b_1 - d_1 a_2
| |
| \,</math>
| |
| | |
| and the effect is that <math>x_1</math> has been eliminated from the second equation. Using a similar tactic with the '''modified''' second equation on the third equation yields:
| |
| | |
| :<math>
| |
| (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
| |
| \,</math>
| |
| | |
| :<math>
| |
| (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.
| |
| \,</math>
| |
| | |
| This time <math>x_2</math> was eliminated. If this procedure is repeated until the <math>n^{th}</math> row; the (modified) <math>n^{th}</math> equation will involve only one unknown, <math>x_n</math>. This may be solved for and then used to solve the <math>(n - 1)^{th}</math> 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:
| |
| | |
| :<math>\tilde a_i = 0\,</math>
| |
| | |
| :<math>\tilde b_1 = b_1\,</math>
| |
| :<math>\tilde b_i = b_i \tilde b_{i - 1} - \tilde c_{i - 1} a_i\,</math>
| |
| | |
| :<math>\tilde c_1 = c_1\,</math>
| |
| :<math>\tilde c_i = c_i \tilde b_{i - 1}\,</math>
| |
| | |
| :<math>\tilde d_1 = d_1\,</math>
| |
| :<math>\tilde d_i = d_i \tilde b_{i - 1} - \tilde d_{i - 1} a_i.\,</math>
| |
| | |
| To further hasten the solution process, <math>\tilde b_i</math> may be divided out (if there's no division by zero risk), the newer modified coefficients notated with a prime will be:
| |
| | |
| :<math>a'_i = 0\,</math>
| |
| | |
| :<math>b'_i = 1\,</math>
| |
| | |
| :<math>c'_1 = \frac{c_1}{b_1}\,</math>
| |
| :<math>c'_i = \frac{c_i}{b_i - c'_{i - 1} a_i}\,</math>
| |
| | |
| :<math>d'_1 = \frac{d_1}{b_1}\,</math>
| |
| :<math>d'_i = \frac{d_i - d'_{i - 1} a_i}{b_i - c'_{i - 1} a_i}.\,</math>
| |
| | |
| This gives the following system with the same unknowns and coefficients defined in terms of the original ones above:
| |
| | |
| :<math>\begin{array}{lcl}
| |
| b'_i x_i + c'_i x_{i + 1} = d'_i \qquad &;& \ i = 1, \ldots, n - 1 \\
| |
| b'_n x_n = d'_n \qquad &;& \ i = n. \\
| |
| \end{array}
| |
| \,</math>
| |
| | |
| 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:
| |
| | |
| :<math>x_n = d'_n/b'_n\,</math>
| |
| | |
| :<math>x_i = (d'_i - c'_i x_{i + 1})/b'_i \qquad ; \ i = n - 1, n - 2, \ldots, 1.</math>
| |
| | |
| ==Variants==
| |
| In some situations, particularly those involving [[periodic boundary conditions]], a slightly perturbed form of the tridiagonal system may need to be solved:
| |
| | |
| :<math>
| |
| \begin{align}
| |
| 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{align}
| |
| </math> | |
| | |
| 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. The method requires solving a modified non-cyclic version of the system for both the input and a sparse corrective vector, and then combining the solutions. This can be done efficiently if both solutions are computed at once, as the forward portion of the pure tridiagonal matrix algorithm can be shared.
| |
| | |
| 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 equation discretized into block tridiagonal|Poisson problem]]). Simplified forms of Gaussian elimination have been developed for these situations{{Citation needed|date=March 2011}}.
| |
| | |
| 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.
| |
| | |
| ==References==
| |
| | |
| <references />
| |
| *{{cite book|author=Conte, S.D., and deBoor, C.|year=1972|title=Elementary Numerical Analysis|publisher= McGraw-Hill, New York|isbn= 0070124469}}
| |
| *{{CFDWiki|name=Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)}}
| |
| *{{Cite book|last1=Press|first1=WH|last2=Teukolsky|first2=SA|last3=Vetterling|first3=WT|last4=Flannery|first4=BP|year=2007|title=Numerical Recipes: The Art of Scientific Computing|edition=3rd|publisher=Cambridge University Press| publication-place=New York|isbn=978-0-521-88068-8|chapter=Section 2.4|chapter-url=http://apps.nrbook.com/empanel/index.html?pg=56|postscript=<!-- Bot inserted parameter. Either remove it; or change its value to "." for the cite to end in a ".", as necessary. -->{{inconsistent citations}}}}
| |
| | |
| ==External links==
| |
| *[http://www.water.tkk.fi/wr/kurssit/Yhd-12.122/www_book/sgh_422b.htm Example with VBA code]
| |
| | |
| {{Numerical linear algebra}}
| |
| | |
| {{DEFAULTSORT:Tridiagonal Matrix Algorithm}}
| |
| [[Category:Numerical linear algebra]]
| |