2.1.7 Sistemas Lineares — Métodos Diretos — O Algoritmo de Thomas

1. O algoritmo de Thomas

Seja um sistema tri-diagonal de \(n \) incógnitas que pode ser escrito como


\( a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i \)

onde \(a_1 = 0 \) e \(c_n = 0 \) . Na forma matricial, este sistema é escrito como

\( \left[
\begin{array}{cccccc}
b_1 & c_1 & 0 & 0 & \cdots & 0 \\
a_2 & b_2 & c_2 & 0 & \cdots & 0 \\
0 & a_3 & b_3 & \ddots & ~ & ~ \\
\vdots & ~ & \ddots & ~ & \ddots & c_{n-1} \\
0 & ~ & ~ & ~ & a_n & b_n
\end{array}
\right] \)

Para a resolução de tal sistema, definimos


\( c_i’ = \left\{ \begin{array}{ll}
\dfrac{c_1}{b_1} & \textrm{para i = 1 }\\
~ \\
\dfrac{c_i}{b_i – c_{i-1}’ a_i} & \textrm{para i = 2, 3, … , n-1}\\
\end{array} \right. \)

e

\( d_i’ = \left\{ \begin{array}{ll}
\dfrac{d_1}{b_1} & \textrm{para i = 1 }\\
~ \\
\dfrac{d_i – d_{i-1}’ a_i}{b_i – c_{i-1}’ a_i} & \textrm{para i = 2, 3, … , n-1 }\\
\end{array} \right. \)

a solução \(x \) deste sistema será

\( x_i = \left\{ \begin{array}{ll}
d_n’ & \textrm{para i = n }\\
~ \\
d_i’ – c_i’ x_{i+1} & \textrm{para i = n-1, n-2, … , 1 }\\
\end{array} \right. \)

2. Implementação

def thomas(A, b):
    """
    Return a list with the solution of linear system.
 
    thomas(A, b)
 
    INPUT:
    * A: a matrix of coefficients
    * b: a list, relative a vector of independent terms
 
    return: a list with unknown vector 'x' of system A*x = b.
 
    Author: Pedro Garcia [sawp@sawp.com.br]
    see: http://www.sawp.com.br
 
    License: Creative Commons
             http://creativecommons.org/licenses/by-nc-nd/2.5/br/
 
    Sep 2010
    """
    n = len(b)
    terms = deepcopy(b)
    x = []
    d = []
    upper = []
    lower = []
 
    for i in xrange(n):
        d.append(A[i][i])
        x.append(0.0)
        upper.append(0.0)
        lower.append(0.0)
    for i in xrange(n - 1):
        upper[i] = A[i][i + 1]
        lower[i] = A[i + 1][i]
 
    for i in xrange(1, n):
        d[i] = d[i] - lower[i - 1] * upper[i - 1] / d[i - 1]
        terms[i] = terms[i] - lower[i - 1] * terms[i - 1] / d[i - 1]
    x[n - 1] = terms[n - 1] / d[n - 1]
 
    # (regressive replacement)
    for i in xrange(n - 2, -1, -1):
        x[i] = (terms[i] - upper[i] * x[i + 1]) / d[i]
    return x

Estas funções estão disponíveis em http://www.sawp.com.br/code/linear/thomas.py, assim como um exemplo de como utilizá-las.

 

3. Copyright

Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.

References

[1] Anthony Ralston and Philip Rabinowitz, A First Course in Numerical Analysis (2nd ed.), McGraw-Hill and Dover, (2001), 423-424.
[2] Lars Elden, Numerical Linear Algebra and Applications in Data Mining (Preliminary version), http://www.mai.liu.se/~ laeld/, (2005)
[3] http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm