Resolução numérica de EDOs de 2ª ordem

Equação diferencial ordinária e seus métodos numéricos de resolução:

Uma equação que contém derivadas de uma função desconhecida é chamada de equação diferencial. Uma equação diferencial com apenas uma variável independente é denominada equação diferencial ordinária (EDO). A ordem de uma EDO é definida pela maior derivada envolvida na equação.

Na solução de uma EDO, procura-se a função que satisfaz a equação analisada. Analiticamente, existem vários métodos de resolução de EDOs. A resolução numérica é utilizada quando temos uma EDO de difícil ou demorada resolução analítica e em problemas de análise, simulação e controle de processos em engenharia, por exemplo.

Numericamente, é possível chegar à solução de uma EDO por diferentes métodos. Há procedimentos que envolvem abordagens de passo simples (a solução y_i+1 no ponto seguinte x_i+1 é obtida a partir da solução y_i conhecida no ponto atual x_i) e outros que utilizam uma abordagem multipasso (a solução em x_i+1 é obtida a partir de soluções conhecidas em vários pontos anteriores). Na solução em cada passo, ainda podem ser usados dois tipos de métodos, explícito e implícito. Métodos explícitos usam fórmulas explícitas para o cálculo da função no próximo valor da variável independente; há apenas variáveis conhecidas no lado direito da equação:

Já nos métodos implícitos, a incógnita aparece de ambos os lados da equação:

De forma geral, métodos implícitos oferecem maior precisão aos resultados que os explícitos, mas exigem maior esforço computacional.

Como ocorre em qualquer procedimento de cálculo numérico, existem erros associados aos métodos de solução numérica de EDOs: erros truncamento e erros de arredondamento. Erros de arredondamento são devido ao procedimento usado por computadores, erros de truncamento ocorrem devido à natureza aproximada dos métodos de solução numérica.

Dois problemas são estudados na resolução de EDOs, os problemas de valor inicial (PVI) e os problemas de valor de contorno (PVC). Um PVI é uma equação diferencial acompanhada do valor da função objetivo em determinado ponto, chamado de valor inicial ou condição inicial. Para a resolução de EDOs de ordens superiores, é necessário conhecer a condição inicial em valores diferentes da variável independente; esses são os PVCs. Suas restrições são geralmente especificadas nos pontos finais do domínio da solução. O número de restrições conhecidas acompanha a ordem da equação (uma equação de n-ésima ordem exige n restrições para a sua resolução).

Para a resolução de PVIs, existem métodos como Euller, Euller modificado, ponto central, Runge-Kutta, multipasso, preditor-corretor. Para os PVCs, os métodos mais conhecidos são método do tiro e método das diferenças finitas. Este tópico tem como foco a explicação do funcionamento e uma aplicação do método do tiro.

Método do tiro:

O método do tiro é um método numérico utilizado para a resolução de PVCs. Nesse método, um PVC é convertido em dois PVIs pela transformação de uma EDO de segunda ordem em doas EDOs de primeira ordem, cada uma com a sua condição inicial. Uma EDO de segunda ordem é convertida em duas EDOs de primeira ordem pela mudança de variável, da seguinte forma:

É introduzida ao problema uma nova variável (w), restando:

Assim, temos:

As condições iniciais do problema são y(a) = Y_a e y(b) = Y_b.

Na solução pelo método do tiro, é realizada uma tentativa para a inclinação (derivada primeira) em x = a, onde temos uma condição inicial Y_a e queremos chegar à segunda condição Y_b.

FIGURA 1: Ilustração do método do tiro.

Fonte: GILAT, Amos. Métodos Numéricos para Engenheiros e Cientistas. Página 411.

Inserida a estimativa, resolve-se as equações (6) e (7) para a obtenção da solução numérica em x = b, y_b1. Esta resolução é feita utilizando-se métodos de resolução de PVIs. Se a solução y_b1 é suficientemente próxima da solução y_b, obedecendo-se uma tolerância, a solução do problema foi obtida. Se não, obtém-se uma nova estimativa até que se chegue a uma solução y_b2 e faz-se novamente a análise em relação a y_b. Quando a tolerância de proximidade com y_b é satisfeita, ou seja, quando y_bi é numericamente igual a y_b, tem-se a solução do problema.

Quando trabalhamos com o MATLAB, podemos implementar uma sub-rotina que faça a análise do erro em relação a y_b e em seguida obter uma estimativa de inclinação (dy/dx) que zere esse erro, chegando assim rapidamente a um valor adequado.

Sugestões:

Para maior aprofundamento teórico nos temas discutidos, recomenda-se a leitura dos capítulos 8 e 9 do livro Métodos Numéricos para Engenheiros e Cientistas (Amos Gilat) e dos capítulos 25 e 27 do livro Métodos Numéricos na Engenharia (Steven Chapra).

Para maior entendimento das funções residentes do MATLAB (ode45) utilizadas para a resolução dos problemas a seguir, recomenda-se a utilização das funções doc ou help do software. Para tal, digite doc ode45 ou help ode45 na janela de comandos e pressione Enter.

Aplicações:

Para exemplificar a utilização do método do tiro em aplicações de engenharia, foram resolvidos dois problemas: a questão da segunda avaliação de ENGD04 (2017.1) que solicita a utilização do método do tiro, cuja solução não se encontra disponível até a data desta postagem e a questão 9.15 do livro-texto Métodos Numéricos para Engenheiros e Cientistas (Amos Gilat). Vejamos:

Questão 2 – 2ª avaliação 2017.1) O escoamento de um fluido sobre um plano inclinado de inclinação θ pode ser representado por uma equação diferencial de segunda ordem dada por:

Onde u(y) é a velocidade do fluido ao longo da lâmina de líquido, ρ a massa específica, g a gravidade µ a viscosidade do fluido. Dado que u(0) = 0 e u’(h) = 0, ρ = 995,7kg/m³, g = 10m/s², µ = 0,798x10³kg/ms, use o método do tiro para encontrar o perfil de u(y) para θ = π/4 e u’(0) = 0,1. Assumir como critério E = û’(h) – u’(h), onde û’(h) é a derivada da velocidade u em relação a y quando y é máximo (ponto h) obtida pelo método do tiro. Considerar h = 0,1m.

Na resolução pelo método do tiro, faremos a mudança de variável, foi explicado acima. Dessa forma, o problema será dividido em dois PVIs, que serão resolvidos pela função residente do MATLAB ode45, que resolve o PVI pelo método de Runge-Kutta de 4ª ordem com erro de truncamento de 5ª ordem. A resolução deste problema foi dividida em três códigos: principal (metodotiro_main), erro (metodotiro_erro) e divisão do PVC em dois PVIs (metodotiro). O programa executável é o principal. É necessário que todos os programas estejam salvos na mesma pasta. Vejamos a resolução:

Programa principal (metodotiro_main):

FIGURA 2: Programa principal

No programa principal, primeiramente são declaradas as condições de contorno, os limites inferior e superior da variável independente (nesse caso o y – altura) e a estimativa inicial do método do tiro, que é a derivada em y = 0 (limite inferior).

Em seguida, é obtida a função F em função da derivada du/dy (que iguala-se a W quando fazemos a mudança de variável) através da sub-rotina erro, que devolve ao programa principal uma função F(W).

Sub-rotina – obtenção da função erro (metodotiro_erro):

FIGURA 3: Sub-rotina para o erro

Neste código é obtida a função erro, que será zerada no código principal pela função residente fsolve, com estimativa inicial W₀. Esta sub-rotina recebe os parâmetros do programa principal, resolve um PVI inicial pela função ode45 utilizando como auxílio a outra sub-rotina, que faz a divisão do PVC em PVIs. Este programa retorna uma função erro em função de W. A análise do erro é realizada em função da derivada de u no último ponto do intervalo de y (u’(h)), obedecendo à condição de contorno dada.

Sub-rotina – divisão do PVC em dois PVIs:

FIGURA 4: Sub-rotina para o mecanismo do método do tiro

Neste programa são definidas as constantes da EDO de segunda ordem, é feita a mudança de variável (mecanismo principal do método do tiro). O programa retorna uma variável dZ/dy, que é um vetor bidimensional, como visto acima.

No retorno ao programa principal, é usada a função fsolve para zerar a função erro usando a estimativa inicial W₀; W é o valor no qual o erro é mínimo. Em seguida, resolve-se os PVIs através da função ode45, que retorna o vetor de variável independente (y) e uma matriz (Z). Essa matriz é composta por várias linhas e duas colunas, sendo a coluna 1 (Z(:,1)) correspondente a u e a coluna 2 (Z(:,2)) correspondente a du/dy (W). Como estamos interessados no comportamento de u em função de y, devemos plotar Z(:,1) vs y. Este procedimento é feito na FIGURA 1. O resultado gráfico obtido é o seguinte: