Algoritmo de Thomas

Origem: Wikipédia, a enciclopédia livre.
Ir para: navegação, pesquisa

O Algorítmo de Thomas (ou Algoritmo de matriz tridiagonal), é um método algébrico que consiste numa simplificação da Eliminação Gaussiana para resolução de sistemas de equações tridiagonais.

Introdução[editar | editar código-fonte]

Uma matriz tridiagonal é uma matriz quadrada onde apenas os elementos da diagonal principal (main) e as que estão acima (up) e abaixo (down) a ela são não nulas. Quando a matriz é tridiagonal, torna-se um desperdício computacional armazenar os zeros. Além disso vemos que eles são desnecessários uma vez que nunca serão utilizados para a solução do sistema. Pensando nisso, Llewellyn Thomas propôs um algoritmo que requer um custo computacional inferior aos métodos de eliminação. Este algoritmo ficou conhecido como Algoritmo de Thomas, o qual requer apenas 8n-7 operações, sendo 3(n-1) operações para a fatorização e 5n-4 operações para o procedimento de substituição.

Método[editar | editar código-fonte]

Uma matriz tridiagonal pode ser escrita na forma:


a_i x_{i - 1}  + b_i x_i  + c_i x_{i + 1}  = d_i , \,\!


onde  a_1  = 0\, y  c_n = 0\, .

e cuja representação na forma matricial se dá por:


\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}
.


O primeiro passo consiste em alterar os coeficientes da seguinte forma:

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


Marcando com o superíndice ' os coeficientes recém alterados.


Da mesma forma faz-se:

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



A solução é então obtida pela substituição de volta:


x_n = d'_n\,
x_i = d'_i - c'_i x_{i + 1} \qquad ; \ i = n - 1, n - 2, \ldots, 1.


Implementações[editar | editar código-fonte]

Implementação em C[editar | editar código-fonte]

A seguinte função em C resolverá o sistema (embora sobreescreverá o vetor de entrada c no processo). Note que o subíndices estão baseados no zero, ou seja, n = 0, 1, \dots, N - 1 onde N é o número de equações:

void solve_tridiagonal_in_place_destructive(float x[], const size_t N, const float a[], const float b[], float c[]) {
size_t n;
 
/*
resolve Ax = v onde A é uma matriz tridiagonal composta pelos veores a, b, c
note que o conteúdo do vetor de entrada c será modificado, tornando esta uma função de um único tempo de uso
x[] - inicialmente contém o vector de entrada v e retorna a solução x. indexados por[0, ..., N - 1]
N - número de equações
a[] - subdiagonal (diagonal abaixo da diagonal principal) -- indexados por [1, ..., N - 1]
b[] - matriz principal, indexados por [0, ..., N - 1]
c[] - superdiagonal (diagonal acima da diagonal principal) -- indexedos por [0, ..., N - 2]
*/
 
c[0] = c[0] / b[0];
x[0] = x[0] / b[0];
 
/* loop de 1 a N - 1 inclusive */
for (n = 1; n < N; n++) {
float m = 1.0f / (b[n] - a[n] * c[n - 1]);
c[n] = c[n] * m;
x[n] = (x[n] - a[n] * x[n - 1]) * m;
}
 
/* loop de N - 2 a 0 inclusive */
for (n = N - 1; n-- > 0; )
x[n] = x[n] - c[n] * x[n + 1];
}



A variante seguinte preserva o sistema de equações para reutilização em outras funções. As chamadas são feitas para a biblioteca para reservar mais espaço. Outras variantes usam um ponteiro para a memória.

void solve_tridiagonal_in_place_reusable(float x[], const size_t N, const float a[], const float b[], const float c[]) {
size_t n;
 
/* alocar espaço scratch */
float * const cprime = malloc(sizeof(float) * N);
 
cprime[0] = c[0] / b[0];
x[0] = x[0] / b[0];
 
/* loop de 1 a N - 1 inclusive */
for (n = 1; n < N; n++) {
float m = 1.0f / (b[n] - a[n] * cprime[n - 1]);
cprime[n] = c[n] * m;
x[n] = (x[n] - a[n] * x[n - 1]) * m;
}
 
/* loop de N - 2 a 0 inclusive */
for (n = N - 1; n-- > 0; )
x[n] = x[n] - cprime[n] * x[n + 1];
 
/* espaço livre scratch  */
free(cprime);
}



Implementação em Python[editar | editar código-fonte]

A seguinte implementação usa linguagem de programação Python. Novamente, os subíndices são basados no zero (i = 0, 1, \dots, n-1 donde n é o número de incógnitas).

def TDMASolve(a, b, c, d):
n = len(d)#n em números é linhas
 
# Modifica o primeiro coeficiente de cada linha
c[0] /= b[0] #Risco de divisão por zero.
d[0] /= b[0]
 
for i in xrange(1, nmax):
ptemp = b[i] - (a[i] * c[i-1])
c[i] /= ptemp
d[i] = (d[i] - a[i] * d[i-1])/ptemp
 
#Substituição de volta
x = [0 for i in xrange(nmax)]
x[-1] = d[-1]
for i in range(-2,-nmax-1,-1):
x[i] = d[i] - c[i] * x[i+1]
 	return x


Implementação em Matlab[editar | editar código-fonte]

Em Matlab o algorítmo fica como segue. Esta vez os vetores estão baseados no um, deste modo i = 1, 2, \dots, n com n que é o número de incógnitas.

function x = TDMAsolver(a,b,c,d)
%a, b, c são os vetores colunas para a compressão da matriz tridiagonal, d é o vetor à direita
% N é o número de linhas
N = length(d);
 
% Modifica o primeiro coeficiente de cada linha
c(1) = c(1) / b(1); % Risco de divisão por zero.
d(1) = d(1) / b(1); 
 
for n = 2:1:N
    temp = b(n) - a(n) * c(n - 1);
    c(n) = c(n) / temp;
    d(n) = (d(n) - a(n) * d(n - 1)) / temp;
end
 
% Agora voltando a substituição.
x(N) = d(N);
for n = (N - 1):-1:1
    x(n) = d(n) - c(n) * x(n + 1);
end
end


Implementação em Fortran 90[editar | editar código-fonte]

Fortran usa também uma nomenclatura baseado no um, ou seja i = 1, 2, \dots, n sendo n o número de incógnitas.

Algumas vezes no é desejável que o programa sobreescreva os coeficientes (por exemplo para resolver diversos sistemas que só difierem no térmo independente), assim que esta implementação mantém ditos coeficientes.

      subroutine solve_tridiag(a,b,c,d,x,n)
      implicit none
!	a - sub-diagonal (diagonal abaixo da diagonal principal)
!	b - diagonal principal
!	c - sup-diagonal (diagonal acima da diagonal principal)
!	d - parte à direita
!	x - resposta
!	n - número de equações
 
        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
 
! inicializar c-primo e d-primo
        cp(1) = c(1)/b(1)
        dp(1) = d(1)/b(1)
! resolver para vetores c-primo e d-primo
         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
! inicializar x
         x(n) = dp(n)
! resolver para x a partir de vetores c-primo e d-primo
        do i = n-1, 1, -1
          x(i) = dp(i)-cp(i)*x(i+1)
        end do
 
    end subroutine solve_tridiag



Exemplo[editar | editar código-fonte]

Seis molas com diferentes constantes k_i e comprimentos L, na condição de repouso são amarradas em série. O ponto final B é então deslocado de tal forma que a distância entre os pontos A e B seja L = 1,2m. Determine as posições x_1, x_2, ..., x_5 dos pontos finais das molas.


As constantes das molas e os seus comprimentos em repouso são:


Mola k (kN/m) L (m)
1 8 0,18
2 9 0,22
3 15 0,26
4 12 0,19
4 10 0,15
6 18 0,30


SOLUÇÃO

A força F em uma mola é dada por:


F = kδ 

onde k é a constante da mola e δ é sua extensão além do comprimento na condição de repouso. Como as molas estão amarradas em série, a força em todas as molas é a mesma. Assim, é possível escrever 5 equações igualando a força em cada duas molas adjacentes. Por exemplo, a condição que a força na primeira mola é igual à força na segunda mola resulta em:


k_1(x_2 - L_1) = k_2 [(x_2 - x_1) - L_2)]


De forma similar, as outras quatro equações são escritas:


k_2[(x_2 - x_1) - L_2] = k_3 [(x_3 - x_2) - L_3]
k_3[(x_3 - x_2) - L_3] = k_4 [(x_4 - x_3) - L_4]
k_4[(x_4 - x_3) - L_4] = k_5 [(x_5 - x_4) - L_5]
k_5[(x_5 - x_4) - L_5] = k_6 [(L - x_5) - L_6]


As cinco equações formas um sistema tridiagonal, cuja forma matricial é:


\begin{bmatrix}
   {k_1 + k_2} & {-k_2} & { 0 } & { 0 } & { 0 } \\
   {-k_2} & {k_2 + k_3} & {-k_3} & { 0  } & { 0  } \\
   { 0  } & {-k_3} & {k_3 + k_4} & {  -k_4 } & { 0 } \\
   { 0  } & { 0  } & {-k_4} & {k_4 + k_5} & {  -k_5 } \\
   { 0 } & { 0 } & {  0 } & {-k_5} & {k_5 + k_6}\\
\end{bmatrix}
\begin{bmatrix}
   {x_1 }  \\
   {x_2 }  \\
   {x_3 }  \\
   {x_4 }   \\
   {x_5 }  \\
\end{bmatrix}
=
\begin{bmatrix}
   {k_1 L_1 - k_2 L_2 }  \\
   {k_2 L_2 - k_3 L_3 }  \\
   {k_3 L_3 - k_4 L_4 }  \\
   {k_4 L_4 - k_5 L_5 }   \\
   {k_5 L_5 + k_6 L - k_6 L_6 }  \\
\end{bmatrix}


O sistema de equações pode ser resolvido com o uso de uma função criada no MATLAB, chamada 'Tridiagonal', exposta abaixo:

function x = Tridiagonal(A,B)
% A função soluciona um sistema tridiaginal de equações lineares [a][x]=[b]
% usando o algorito de Thomas.
% Variável de entrada:
% A Matriz  de coeficientes.
% B Vetor coluna de constantes.
% Variável de saída.
% x Vetor coluna com a solução.
 
%PASSO 1: 
%Define o vetor d com os elementos da diagonal.
[nR, nC] = size(A);
for i  = 1:nR
   d(i) = A(i,i);
end
 
%Define o vetor ad com os elementos acima da diagonal.
for i  = 1:nR - 1
   ad(i) = A(i,i+1);
end 
 
%Define o vetor bd com os elementos abaixo da diagonal.
for i  = 2:nR
   bd(i) = A(i,i-1);
end 
 
%PASSO 2:
ad(1) = ad(1)/d(1);
B(1) = B(1)/d(1);
 
%PASSO 3:
for i = 2:nR - 1
   ad(i) = ad(i)/(d(i) - bd(i)*ad(i-1));
   B(i) = (B(i) - bd(i)*B(i-1))/(d(i) - bd(i)*ad(i-1));
end
 
%PASSO 4:
B(nR) = (B(nR) - bd(nR)*B(nR-1))/(d(nR) - bd(nR)*ad(nR-1));
x(nR,1) = B(nR);
 
%PASSO 5:
for i = nR-1:-1:1
   x(i,1) = B(i) - ad(i)*x(i+1);
end


Utilizando a função Tridiagonal acima para resolver o sistema de equações do problema:

k_1 = 8000; k_2 = 9000; k_3 = 15000; k_4 = 12000; k_5 = 10000; k_6 = 18000;
L = 1.5; L_1 = 0.18; L_2 = 0.22; L_3 = 0.26; L_4 = 0.19; L_5 = 0.15; L_6 = 0.30;   
a = [k_1 + k_2,-k_2, 0, 0, 0; -k_2, k_2 + k_3, -k_3, 0, 0; 0, -k_3, k_3 + k_4, -k_4, 0; 0, 0, -k_4, k_4 + k_5, -k_5; 
0, 0, 0, -k_5, k_5 + k_6];
b = [k_1 L_1 - k_2 L_2; k_2 L_2 - k_3 L_3; k_3 L_3 - k_4 L_4; k_4 L_4 - k_5 L_5; k_5 L_5 + k_6 L - k_6 L_6];
Xs = Tridiagonal(a,b)


A solução encontrada, quando o arquivo é executado é:


\begin{bmatrix}
   {0.2262}\\
   {0.4872}\\
   {0.7718}\\
   {0.9926}\\
   {1.1795}\\
\end{bmatrix}


Derivação[editar | editar código-fonte]

Algorítmo[editar | editar código-fonte]

O algorítmo pode ser obtido a partir da Eliminação Gaussiana na sua forma genérica. Supondo como incógnitas  x_1,\ldots, x_n e com as equações a serem resolvidas:

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


Modificando a segunda equação a partir da primera obtemos:


(\mbox{equação 2}) \cdot b_1 - (\mbox{equação 1}) \cdot a_2


resultando:


(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
\,

(b_2 b_1 - c_1 a_2) x_2  + c_2 b_1 x_3 = d_2 b_1 - d_1 a_2
\,


e se eliminou x_1 da segunda equação. Se repetimos o processo usando a segunda equação na terceira obtemos:


(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
\,

(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.
\,


Esta vez eliminou-se x_2. O procedimento pode ser repetido até a enésima fila, deixando cada equação com somente duas incógnitas, exceto a última que só terá uma incógnita. Então é imediata a resolução desta e, por consequência, das equações interiores.

Fórmulas para coeficientes[editar | editar código-fonte]

A determinação dos coeficientes nas equações gerais é mais complicada. Examinando el procedimento, se podem definir de forma recursiva:

\tilde a_i = 0\,
\tilde b_1 = b_1\,
\tilde b_i = b_i \tilde b_{i - 1} - \tilde c_{i - 1} a_i\,
\tilde c_1 = c_1\,
\tilde c_i = c_i \tilde b_{i - 1}\,
\tilde d_1 = d_1\,
\tilde d_i = d_i \tilde b_{i - 1} - \tilde d_{i - 1} a_i.\,


Para acelerar la resolução \tilde b_i pode ser dividido (se não houver problemas de divisão por zero e os novos coeficiente, indicados com y los nuevos coeficientes, indicados com aspas, serão:

a'_i = 0\,
b'_i = 1\,
c'_1 = \frac{c_1}{b_1}\,
c'_i = \frac{c_i}{b_i - c'_{i - 1} a_i}\,
d'_1 = \frac{d_1}{b_1}\,
d'_i = \frac{d_i - d'_{i - 1} a_i}{b_i - c'_{i - 1} a_i}.\,


Resultando no seguinte sistema com as mesmas incógnitass e os coeficientes definidos em função dos originais:

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


Onde a última equação só afeta uma incógnita. Resolvendo-a, pode-se resolver a anterior e assim sucessivamente:

x_n = d'_n/b'_n\,
x_i = (d'_i - c'_i x_{i + 1})/b'_i \qquad ; \ i = n - 1, n - 2, \ldots, 1.


Variações[editar | editar código-fonte]

Em alguma situações, em particular aquelas envolvendo condições de fronteira, uma ligeiramente perturbada forma do sistema tridiagonal necessita ser resolvido:


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


Neste caso, usa-se a fórmula de Sherman-Morrison para se evitar as operações adicionais da Eliminação Gaussiana e ainda usar o Algoritmo de Thomas. Assim, resolve-se um problema na forma:

 
(A'+uv^T)x = d

onde


u^T: [-b_1 0 0 ... 0 -a_1 \ b_1].


Considere um sistema tridiagonal ligeiramente perturbado A' da mesma forma do exemplificado acima. A solução para esse sistema é obtida resolvendo



A'y = d, A'q = u


E escreva x como:


x = y-[(v^T y)/(1+(v^Tq))]q


Em outras situações, o sistema de equações pode ser tridiagonal de bloco, com submatrizes menores, dispostas como os elementos individuais do sistema de matriz acima. Formas simplificadas de eliminação de Gauss foram desenvolvidas para estas situações.

O livro didático de Matemática Numérica por Quarteroni, Sacco e Saleri, apresenta uma versão modificada do algoritmo que evita algumas das divisões (fazendo uso de multiplicações), o que é benéfico em algumas arquiteturas de computadores.


Referências Bibliográficas[editar | editar código-fonte]

  • Thomas, L.H. Elliptic Problems in Linear Differential Equations over a Network. Columbia University, New York. 1949.
  • Conte, S.D., and deBoor, C. Elementary Numerical Analysis. McGraw-Hill (ed.), New York. 1972. ISBN = 0070124469.
  • Este artigo inclui texto do artigo Tridiagonal matrix algorithm - TDMA (Thomas algorithm) publicado com licência GNU em CFD online wiki
  • Press,W.H., Teukolsky,S.A., Vetterling,W.T., Flannery,B.P. Numerical Recipes: The Art of Scientific Computing. 3rd ed, Cambridge University Press, New York, 2007. Section 2.4. ISBN = 978-0-521-88068-8
  • Algorítmos de Resolução de Sistemas Lineares com Estrutura Tridiagonal. Instituto de Ciências Exatas - ICE- Universidade Federal de Itajubá (UNIFEI) e Instituto de Engenharia de Sistemas e Tecnologias da Informação - IESTI - Universidade Federal de Itajubá (UNIFEI). 2010.
  • Gilat, A. and Subramaniam, V.. Métodos Numéricos para Engenheriros e Cientistas: Uma introdução com aplicação usando o MATLAB. [S.l.]: Bookman (ed.), 2000.