Algoritmo de Thomas
# 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:
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.

### 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).

### 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.

### 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.

### Implementação em Javascript
## 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:
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:
Utilizando a função Tridiagonal acima para resolver o sistema de equações do problema:
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 !CS1 manut: Nomes múltiplos: lista de autores (link)