2.1.7 Sistemas Lineares — Métodos Diretos — O Algoritmo de Thomas
1. O algoritmo de Thomas
Seja um sistema tri-diagonal de incógnitas que pode ser escrito como
$$ a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i $$
onde e . 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 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 </p>
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/.