Métodos das Diferenças Finitas

Teoria Completa

Introdução

E aí meus jovens, tudo bem? Bom, desde o início do estudo de Resoluções Numéricas de EDO, só que o que vimos até agora foram métodos que nós usávamos para resolver problemas de valores iniciais, com essa carinha aqui: y ' ' = f x , y , y ' ' y a = α y ' a = β

Basicamente, temos condições de contorno definidas no mesmo ponto, nesse caso, o ponto x = a .

A partir de agora, a conversa muda e o negócio vai ficar animado!

Vamos ver agora um método para resolver problemas desse tipo aqui:

y ' ' = f x , y , y ' ' y a = α y b = β

Nesse caso, as nossas condições de contorno são definidas pontos diferentes, a e b . Chamamos esse tipo de problema de Problema de Valor de Contorno, ou PVC.

Método das Diferenças Finitas

Tipo assim, o que queremos com esse método é substituir cada derivada por uma aproximação daquelas que vimos lá no começo desse assunto onde aprendemos a derivar numericamente!

Então se você não ta lembrando muito bem recomendo fortemente que dê uma olhada lá!

Podemos começas a entender esse problema pelo teorema de Taylor, podemos escrever o desenvolvimento em série de uma função y ( x ) em torno do ponto x   =   x i como:

y x = y x i + x - x i y ' x i + ( x - x i 2 ) 2 ! y ' ' x i + …  

E o que acontece depois disso? Vamos fazer umas substituições encima dessa série, mas antes de começas a parte técnica eu vou te contar o objetivo de fazer tudo isso.

Tá ligado que nós temos um objetivo aqui, que é fugir do problema de resolver aquelas “EDO”zonas para um problema de resolver um sistema algébrico.

Com esse método dividimos o intervalo a , b cujos extremos são os pontos das condições de contorno, em n partes iguais de tamanho h . Ou seja:

h = b - a n ⇔ n = b - a h

No final das contas o método vai nos dar uma aproximação do valor de y x para cada x i em que dividimos o intervalo a , b .

Fazendo algumas modificações na série e substituindo x por x i + 1 e brincando um pouco com os números encontramos isso aqui quando queremos uma Diferença finita ascendente, que mostra a derivada primeira de uma função aproximada por:

y i ' ≅ y ( x i + 1 ) - y ( x i ) h

Agora se mudamos x por x i - 1 temos a chamada Diferença finita descendente:

y i ' ≅ y ( x i ) - y ( x i - 1 ) h

Chamamos essas duas de aproximação de primeira ordem, onde o limitante do erro é dado por:

E ≤ h 2 ∙ max x ∈ a , b ⁡ y ' ' x

Uma terceira aproximação para y ’ ( x ) pode ser obtida subtraindo-se a expressão para   y ( x i - 1 ) da expressão para y ( x i + 1 ) , ou seja:

Diferença Central:

y i ' ≅ y ( x i + 1 ) - y ( x i - 1 ) 2 h

Essa é uma aproximação de segunda ordem, onde o limitante do erro é dado por

E ≤ h 2 6 ∙ max x ∈ a , b ⁡ y 3 x

Já para a derivada segunda y ' ' x usaremos apenas a fórmula da Diferença Central para a Derivada Segunda:

y i ' ' ≅ y i + 1 - 2 y i + y i - 1 h 2

Ela também é uma aproximação de segunda ordem e o limitante do erro é dado por:

E ≤ h 2 12 ∙ max x ∈ a , b ⁡ f 4 x

Obs.: Como você deve imagina, quanto maior o grau de aproximação melhor é o resultado.

Obs².: Se a EDO for linear, o método vai torna-la um sistema linear que pode ser resolvido por Eliminação de Gauss, por exemplo. Se a EDO for não linear, resultará em um sistema não linear que pode ser resolvido por Newton. Nesse tópico nos focaremos apenas nas lineares.

Exemplo

Vamos ver direto a forma geral. Resolve aí:

y ' ' = p x y ' + q x y + r x y a = α y b = β

Primeiro dividiríamos o intervalo a , b em n intervalos iguais de comprimento h . Daí, teríamos:

x i = a + h i

Agora, usando os métodos de diferença central, substituímos a EDO:

y ' ' = p x y ' + q x y + r x

y i + 1 - 2 y i + y i - 1 h 2 = p x i y i + 1 - y i - 1 2 h + q x i y i + r x i

Manipulando isso aí e substituindo os x i , chegamos em um sistema A y = b tal que:

A = 2 + h 2 q x 1 - 1 + h 2 p x 1 0 ⋯ 0 - 1 - h 2 p x 2 2 + h 2 q x 2 - 1 + h 2 p x 2 ⋯ ⋮ 0 - 1 - h 2 p x 3 ⋱ ⋱ 0 ⋮ ⋮ ⋱ 2 + h 2 q x n - 2 - 1 + h 2 p x n - 2 0 ⋯ 0 - 1 - h 2 p x n - 1 2 + h 2 q x n - 1

y = y 1 y 2 ⋮ y n - 2 y n - 1                                                   b = - h 2 r x 1 + 1 + h 2 p x 1 y 0 - h 2 r x 2 ⋮ - h 2 r x n - 2 - h 2 r x n - 1 + 1 - h 2 p x n - 1 y n

Bom, a menos que você tenha uma memória sobrenatural não recomendo tentar gravar isso! Essa forma geralzona vai ser útil pra escrever nosso algoritmo. Quando estivermos resolvendo o problema na mão, vamos montar esse sistema na hora. Vai ficar mais fácil de entender nos exercícios!

Algoritmo

Se prepara que esse é grande em! Queremos resolver um sistema do tipo:

y ' ' = p x y ' + q x y + r x

a ≤ x ≤ b

y a = α

y b = β

Passo 1. Dados iniciais (o que seu programa vai receber de você):

  1. Número de subintervalos n
  2. Intervalo a , b
  3. Condições de contorno α e β

Passo 2. h = b - a n

Passo 3. Defina as variáveis abaixo.

x = a + h

a 1 = 2 + h 2 q x

b 1 = - 1 + h 2 p x

d 1 = - h 2 r x + 1 + h 2 p x α

Passo 4. Para i variando de 2 até n - 2 repita o passo 5.

Passo 5. Defina as variáveis abaixo.

x = a + h i

a i = 2 + h 2 q x

b i = - 1 + h 2 p x

c i = - 1 - h 2 p x

d i = - h 2 r x

Passo 6. Defina as variáveis abaixo.

x = b - h

a n - 1 = 2 + h 2 q x

c n - 1 = - 1 - h 2 p x

d n - 1 = - h 2 r x + 1 - h 2 p x β

Passo 7. Defina as variáveis abaixo.

l 1 = a 1

u 1 = b 1 a 1

z 1 = d 1 l 1

Passo 8. Para i variando de 2 até n - 2 repita o passo 9.

Passo 9. Defina as variáveis abaixo.

l i = a i - c i u i - 1

u i = b i l i

z i = d i - c i z i - 1 l i

Passo 10. Defina as variáveis abaixo.

l n - 1 = a n - 1 - c n - 1 u n - 2

z n - 1 = d n - 1 - c n - 1 z n - 2 l n - 1

Passo 11. Defina as variáveis abaixo.

y 0 = α

y n = β

y n - 1 = z n - 1

Passo 12. Para i variando de n - 2 até 1 repita o passo 13.

Passo 13. Defina a variáveis abaixo.

x i = z i - u i y i + 1

Passo 14. Retorne os y i .

E é isso aí. O passo 7 até o 13 são uma resolução mais eficiente para o sistema. Você também pode resolver usando qualquer outro método que aprendemos!

Pronto galera, agora vamos praticar nos exercícios!

Expandir

Exercício Resolvido #1

Márcia A. Gomes Ruggiero e Vera Lúcia da Rocha Lopes, Cálculo Numérico – Aspectos Teóricos e Computacionais, 2ª ed. São Paulo: MAKRON Books, 1996, pag 374, exercício 24

Resolva pelo método das diferenças finitas, o PVC:

y ' ' + 2 y ' + y = x y 0 = 2 y 1 = 0

usando h = 0.25

MOSTRAR SOLUÇÃO COMPLETA

Passo 1

Temos o seguinte:

y ' ' + 2 y ' + y = x

E teremos n intervalos de modo que:

n = b - a h = 1 0.25 = 4

Então vai ficar assim:

Beleza, agora vamos substituir as derivadas na nossa EDO:

y ' ' + 2 y ' + y = x

y i + 1 - 2 y i + y i - 1 h 2 + 2 y i + 1 - y i - 1 2 h + y i = x i

Multiplicando os dois lados por h 2 :

y i + 1 - 2 y i + y i - 1 + h y i + 1 - y i - 1 + h 2 y i = h 2 x i

Juntando os termos iguais:

y i + 1 1 + h + y i - 2