Eliminação gaussiana e retro-substituição

Em diversos problemas da fisica, como na interpolação por Spline cúbicos, resolução de equações diferenciais, etc, envolvem equações do tipo

$b_{i}^{-}\phi _{i-1}+b_{i}^{0}\phi _{i}+b_{i}^{+}\phi _{i+1}=f_{i}.\;\;\;\;\;\;(1)$

ou na forma matricial

$\mathbf {b\phi } =\mathbf {f} \;.$

A matriz $\mathbf {b}$ e o vetor $\mathbf {f}$ são dados do problema e $\mathbf {\phi }$ representa a solução procurada. A primeira vista, trata-se de um problema simples, bastanto calcular $\mathbf {b} ^{-1}$ para se obter:

$\mathbf {\phi } =\mathbf {b^{-1}f} \;.$

Contudo, alguns aspectos devem ser levados em consideração. Em primeiro lugar, do ponto de vista numérico, a inversão de matrizes é um problema sutil e muito oneroso do ponto de vista computacional (veja Numerical Recipes, por exemplo). Por isto, é fortemente recomendável que se procure algoritmos numéricos bem estabelecidos para a realização da tarefa. Tendo em vista que tais algoritmos estão implementados e testados há várias décadas, isto leva a crer que não há dificuldades na resolução do problema. Entretanto, o custo em tempo é algo não desprezível e pode comprometer seriamente um estudo no qual esta equação tenha que ser resolvida muitas vezes. Por exemplo, se os coeficientes $\;b_{i}^{-}$ , $\;b_{i}^{0}$ e $\;b_{i}^{+}$ , e/ou $\;f_{i}$ variam com o tempo, a matriz $\mathbf {b}$ deve ser invertida a cada passo da evolução.

Entretanto, nos casos em que o sistema de equações envolve apenas os termos em $\;\phi _{i-1}$ , $\;\phi _{i}$ e $\;\phi _{i+1}$ , a matriz $\mathbf {b}$ é tri-diagonal, isto é, apenas os elementos de suas 3 diagonais principais são não-nulos. Neste caso, a soulção do problema é bastante segura, simples e eficiente do ponto de vista numérico.

De fato, o sistema de equações acima sugere que a solução pode ser escrita da seguinte forma:

$\;\phi _{i+1}=\alpha _{i}\phi _{i}+\beta _{i}\;\;\;\;\;\;(2)$

onde $\;\alpha _{i}$ e $\;\beta _{i}$ devem ser determinados a partir do sistema de equações. Substituindo a expressão acima na Eq. (1), temos:

$b_{i}^{-}\phi _{i-1}+(b_{i+1}^{+}\alpha _{i}+b_{i}^{0})\phi _{i}=f_{i}-b_{i}^{+}\beta _{i}$

o que leva a

$\phi _{i}=\gamma _{i}b_{i}^{-}\phi _{i-1}+\gamma _{i}(b_{i}^{+}\beta _{i}-f_{i})\;,\;\;\;\;\;\;(3)$

onde

$\gamma _{i}\equiv -{\frac {1}{b_{i}^{+}\alpha _{i}+b_{i}^{0}}}\;.\;\;\;\;\;\;(4)$

Reescrevendo a Eq. (2) como

$\;\phi _{i}=\alpha _{i-1}\phi _{i-1}+\beta _{i-1}$

podemos comparar esta expressão com a Eq. (3), o que nos leva, finalmente, a

$\;\alpha _{i-1}=\gamma _{i}b_{i}^{-}\;\;\;\;\;\;(5)$

e

$\;\beta _{i-1}=\gamma _{i}(b_{i}^{+}\beta _{i}-f_{i})\;.\;\;\;\;\;\;(6)$

As equações (2) e (4-6) mostram que a solução do problema é bastante simples. Dado o valor de $\;\phi _{N}=q$ , temos $\;\alpha _{N-1}=0$ e $\;\beta _{N-1}=q$ . Uma vez obtidas estas quantidades, calculamos $\;\alpha _{i}$ e $\;\beta _{i}$ , por meio das equações (4-6), de $\;i=N-1$ até $\;i=2$ . Após isto, $\;\phi _{i+1}$ pode ser obtido por meio da equação (2), partindo de $\;i=1$ até $\;i=N-2$ .

É importante notar que a solução do problema só é possível se forem fornecidos os valores de $\;\phi _{1}$ e $\;\phi _{N}$ . A especificação destes valores está associada às condições de contorno, que são determinadas pela física do problema. Na realidade, também é comum se encontrar problemas onde os valores das derivadas são fixados em uma ou ambas das bordas. Se $\;\phi '_{1}=\sigma _{1}$ ou $\;\phi '_{N}=\sigma _{2}$ forem impostos, temos: