Exibir código-fonte para Eliminação gaussiana e retro-substituição

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

<math>
b^-_i\phi_{i-1} + b^0_i\phi_i + b^+_i\phi_{i+1} = f_i.\;\;\;\;\;\; (1)
</math>

ou na forma matricial

<math>
\mathbf{b\phi}=\mathbf{f}\;.
</math>

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

<math>
\mathbf{\phi}=\mathbf{b^{-1}f}\;.
</math>

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 [http://www.nr.com 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 <math>\;b^-_i</math>, <math>\;b^0_i</math> e <math>\;b^+_i</math>, e/ou <math>\;f_i</math> variam com o tempo, a matriz <math>\mathbf{b}</math> deve ser invertida a cada passo da evolução.


Entretanto, nos casos em que o sistema de equações envolve apenas os termos em <math>\;\phi_{i-1}</math>, <math>\;\phi_i</math> e <math>\;\phi_{i+1}</math>, a matriz <math>\mathbf{b}</math> é 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:

<math>
\;\phi_{i+1}=\alpha_i\phi_i+\beta_i\;\;\;\;\;\; (2)
</math>

onde <math>\;\alpha_i</math> e <math>\;\beta_i</math> devem ser determinados a partir do sistema de equações.
Substituindo a expressão acima na Eq. (1), temos:

<math>
b^-_i\phi_{i-1}+(b^+_{i+1}\alpha_i+b^0_i)\phi_i=f_i-b^+_i\beta_i
</math>

o que leva a

<math>
\phi_i=\gamma_i b^-_i\phi_{i-1}+\gamma_i(b^+_i\beta_i-f_i)\;,\;\;\;\;\;\; (3)
</math>

onde

<math>
\gamma_i\equiv-\frac{1}{b^+_i\alpha_i+b^0_i}\;.\;\;\;\;\;\; (4)
</math>

Reescrevendo a Eq. (2) como

<math>
\;\phi_i=\alpha_{i-1}\phi_{i-1}+\beta_{i-1}
</math>

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

<math>
\;\alpha_{i-1}=\gamma_i b^-_i\;\;\;\;\;\; (5)
</math>

e

<math>
\;\beta_{i-1}=\gamma_i(b^+_i\beta_i-f_i)\;.\;\;\;\;\;\; (6)
</math>


As equações (2) e (4-6) mostram que a solução do problema é bastante simples. Dado o valor de <math>\;\phi_N=q</math>, temos <math>\;\alpha_{N-1}=0</math> e <math>\;\beta_{N-1}=q</math>. Uma vez obtidas estas quantidades, calculamos <math>\;\alpha_i</math> e <math>\;\beta_i</math>, por meio das equações (4-6), de
<math>\;i=N-1</math> até <math>\;i=2</math>. Após isto, <math>\;\phi_{i+1}</math> pode ser obtido por meio da equação (2), partindo de <math>\;i=1</math> até <math>\;i=N-2</math>.


É importante notar que a solução do problema só é possível se forem fornecidos os valores de <math>\;\phi_1</math> e <math>\;\phi_N</math>. 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 <math>\;\phi'_1=\sigma_1</math> ou <math>\;\phi'_N=\sigma_2</math> forem impostos, temos:

<math>
\phi'_1=\frac{\phi_2-\phi_1}{X_2-X_1}=\sigma_1
</math>

e

<math>
\phi'_N=\frac{\phi_N-\phi_{N-1}}{X_N-X_{N-1}}=\sigma_2
</math>

o que leva a

<math>
\;\phi_2=\phi_1+\sigma_1(X_2-X_1)
</math>

e

<math>
\phi_N=\phi_{N-1}+\sigma_2(X_N-X_{N-1})\;.
</math>

Usando a equação (2), encontramos:

<math>
\;\alpha_{N-1}=1\;\;\; \mbox{ e }\;\;\; \beta_{N-1}=\sigma_2(X_N-X_{N-1})
</math>

e

<math>
\phi_1=\frac{\beta_1-\sigma_1(X_2-X_1)}{1-\alpha_1}.
</math>

A partir destes resultados, as relações de recorrência podem ser aplicadas.


----

Voltar para o índice de [[Métodos computacionais]].