Algoritmo de Thomas

Em álgebra linear, o Algoritmo de Thomas (ou Algoritmo de matriz tridiagonal), é um método algébrico oriundo de uma simplificação da eliminação gaussiana para resolução de sistemas de equações tridiagonais.

Uma matriz tridiagonal é uma matriz quadrada onde apenas os elementos da diagonal principal e as que estão acima e abaixo a ela são não nulas. Quando a matriz é tridiagonal, torna-se um desperdício computacional armazenar os zeros, já que eles 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

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\,$ e $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

Implementação em C

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

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

Em Matlab o algoritmo 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

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

Implementação em Javascript

/*
 * Exemplo
 * var t = new Thomas ({
 *		a:[],
 *		b:[],
 *		c:[],
 *		d:[],
 * });
 * console.log(t.x);
 */
function Thomas(opt) {
	opt.c[0] /= opt.b[0];
	opt.d[0] /= opt.b[0];
	x=[];
	
	for (n=1; n < opt.b.length - 1 ; n++) {
		temp = opt.b[n] - (opt.a[n] * opt.c[n - 1]);
		opt.c[n] /= temp;
		opt.d[n] = (opt.d[n] - opt.a[n] * opt.d[n - 1]) / temp;
	}

	x[opt.b.length - 1] = opt.d[opt.b.length - 1];
	for (n = opt.b.length - 2; n >= 0; n-- ) {
		x[n] = opt.d[n] - opt.c[n] * x[n + 1];
	}
	return {
		x: x
	};
}

Exemplo

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_{1}-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

Algoritmo

O algoritmo 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{aligned}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{aligned}}

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

A determinação dos coeficientes nas equações gerais é mais complicada. Examinando o 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 a resolução, ${\tilde {b}}_{i}$ pode ser dividido (se não houver problemas de divisão por zero) e os novos coeficientes (indicados com ') 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ógnitas e os coeficientes definidos em função dos originais:

{\begin{array}{lcl}x_{i}+c'_{i}x_{i+1}=d'_{i}\qquad &;&\ i=1,\ldots ,n-1\\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}\,

x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1.

Variações

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

{\begin{aligned}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{aligned}}

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}00...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^{T}q))]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

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. (2000). Métodos Numéricos para Engenheriros e Cientistas: Uma introdução com aplicação usando o MATLAB. [S.l.]: Bookman