Resoluçao numérica de EDOs

Prévia do material em texto

ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Métodos Matemáticos e Computacionais na Engenharia
ENG-D04 (Parte III)
Resolução Numérica de EDOs
DEPARTAMENTO DE ENGENHARIA QUÍMICA ESCOLA POLITÉCNICA
UFBA
Prof. Leizer Schnitman
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
CONTEÚDO PROGRAMÁTICO
Previsão = 10h (5 aulas)
3. Resolução Numérica de EDOs
3.1 Introdução
3.2 Métodos de passo simples
3.2.1 Método de Euler explícito
3.2.2 Método de Euler modificado
3.2.3 Método do Ponto Central
3.2.4 Método de Euler implícito
3.2.5 Método de Range-Kutta
3.3 Métodos de passo múltiplo
3.3.1 Método de Adams-Bathfort
3.3.2 Método de Adams-Moulton
3.3.3 Método de predição-correção
3.4 Recursos do Matlab
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.1 Introdução
Seja uma função y=f(x) expressa na forma:
Sabe-se que: 
Nosso primeiro objetivo é fazer o caminho inverso, ou seja, conhecida a derivada dy/dx, encontramos a função y=f(x)
Assim, separando dy e dx e integrando os dois lados da equação, ficamos com:
A expressão acima resgata a função original, exceto pelo valor da constante 1 que foi perdida, representada agora por uma constante C, chamada de constante de integração.
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Existem, porém, infinitos valores para a constante C. As figuras seguintes ilustram a situação.
A conclusão é que, para se caracterizar a solução de uma equação diferencial, são necessárias algumas condições auxiliares
Por exemplo, pode-se definir a condição inicial. No caso anterior, se fosse definida uma condição inicial x=0 e y=1 (ou f(0)=1), os valores poderiam ser substituídos na equação para determinar o valor de C=1.
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
A solução numérica de uma EDO é formada por um conjunto de pontos discretos que representam a função f(x) de uma forma aproximada
A solução pode ser caracterizada em
 Passos simples: quando a solução no ponto seguinte (xi+1) é calculada com base na solução obtida para o ponto atual (xi)
 Passos múltiplos: quando a solução no ponto seguinte (xi+1) é calculada com base na solução obtida para vários pontos anteriores
 Métodos explícitos: usam uma forma explícita para o cálculo do próximo valor da função, ou seja yi+1 = f(xi, xi+1, yi)
 Métodos implícitos: usam uma forma implícita para o cálculo do próximo valor da função, ou seja yi+1 = f(xi, xi+1, yi+1)
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Quanto ao erro numérico cometido:
 Erro de arredondamento: causado pelo procedimento usado pelo computador
 Erro de truncamento local: associa erro causado pelos métodos numéricos de passo simples
 Erro de truncamento acumulado: é causado pelo acumulo de erros de truncamentos locais gerados nos passos anteriores
 Erro de truncamento global (total): é a soma dos erros de arredondamento e de truncamento
Genericamente, num método explícito de passo simples, dada a solução conhecida no ponto (xi,yi), a solução é calculada usando:
xi+1 = xi + h e yi+1 = yi + Inclinação.h
onde h é a largura do passo de integração e Inclinação é uma constante que estima o valor de dy/dx.
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.2 Métodos de passo simples
3.2.1 Método explícito de Euler (também chamado de Euler progressivo)
A solução é calculada usando: xi+1 = xi + h e yi+1 = yi + Inclinação.h
onde a inclinação é dada pela própria equação diferencial:
Observar que neste contexto a f(x,y) é a própria função da derivada
Note que o método assume que, em uma pequena distância h na vizinhança de (xi,yi), a função tem uma inclinação constante e igual a inclinação neste ponto (xi,yi). Assim:
		xi+1 = xi + h e yi+1 = yi + f(xi,yi).h
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
OBSERVAÇÃO:
Escreva a equação na forma integral:
Assim: 
No módulo anterior, vimos vários métodos de integração numérica e um dos mais simples era o método do retângulo, onde o integrando é representado de forma aproximada por uma constante f(xi,yi).
Com esta abordagem:
Que é a mesma expressão do slide anterior !!
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Exemplo:
Seja equação diferencial de x=0 a x=2,5 e a condição inicial x=0 e y=3
Resolva a EDO fazendo h=0,5; h=0,1 e h=0,01.
Compare com a solução analítica
Ver exemplo Matlab – Euler explícito
Observar que os cálculos serão:
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Análise do erro do método de Euler explícito
Desenvolver !!
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.2.2 Método de Euler modificado
Ë uma técnica explícita de passo simples.
No método de Euler explícito, assumimos que a derivada entre os pontos (xi,yi) e (xi+1, yi+1) era constante e igual a derivada de y(x) no ponto (xi,yi)
Método explícito de Euler
Erro cometido no 1º passo
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
No método de Euler modificado, a inclinação utilizada é modificada buscando incluir a sua variação ao longo do intervalo.
A inclinação utilizada neste método é a média entre a inclinação no início do intervalo e a inclinação no final do intervalo.
A inclinação no início do intervalo e dada por:
Para o cálculo da inclinação no final do intervalo :
 Calcula-se o valor aproximado para yi+1 usando método explícito de Euler
 Estima-se a inclinação do final do intervalo com a substituição do ponto (xi+1,yi+1) na Eq (1)
 Com base no valor médio das inclinações, estima-se yi+1:
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
As figuras seguintes ilustram o esquema do método de Euler modificado:
Algoritmo:
Calcule f(xi,yi)
Calcule o próximo valor da variável independente: xi+1 = xi + h
Estime yi+1 usando o método de Euler explícito: yi+1 = yi + f(xi,yi)h
Calcule f(xi+1,yi+1)
Calcule a inclinação média:
 
Calcule a solução numérica:
Ver exemplo Matlab – Euler modificado
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.2.3 Método do ponto central
Ë uma técnica explícita de passo simples.
No método de Euler explícito, assumimos que a derivada entre os pontos (xi,yi) e (xi+1, yi+1) era constante e igual a derivada de y(x) no ponto (xi,yi).
No método de Euler modificado assumimos que a derivada (inclinação) era dada pela média da inclinação nos pontos (xi,yi) e (xi+1, yi+1).
Neste caso, a inclinação usada é a estimativa da inclinação no ponto central do intervalo, ou seja xm = xi + h/2
Algoritmo:
Calcule f(xi,yi) 
Calcule o ponto médio: xm = xi + h/2
Estime ym usando o método de Euler explícito: ym = yi + f(xi,yi)h/2
Calcule f(xm,ym)
Calcule a solução numérica:
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
As imagens seguintes ilustramo esquema proposto
Calcule f(xi,yi) 
Calcule o ponto médio: xm = xi + h/2
Estime ym usando o método de Euler explícito: ym = yi + f(xi,yi)h/2
Calcule f(xm,ym)
Calcule a solução numérica:
Ver exemplo Matlab – Ponto Central
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Neste caso, porém, a inclinação, ao invés de ser dada pela inclinação no início do intervalo, ela é dada pela inclinação do ponto final do intervalo. A imagem ao lado ilustra o caso.
Com base nisso, o ponto seguinte da equação numérica é calculado por:
xi+1 = xi + h e yi+1 = yi + f(xi+1, yi+1).h
Observar que o termo yi+1 aparece nos dois lados da equação
3.2.4 Método Implícito de Euler
Como todos os métodos de passo simples vistos anteriormente, a solução é calculada como: 
xi+1 = xi + h e yi+1 = yi + Inclinação.h
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Exemplo:
Um componente químico com concentração inicial n0=2000 em t=0, apresenta um decaimento temporal em sua concentração quando exposto ao ar. Esse decaimento ocorre em uma taxa proporcional a sua concentração elevada a 3/2. Simultaneamente, o componente é produzido por outro processo. A equação diferencial que descreve a concentração instantânea do componente é dada por:
Resolva a equação diferencial para obter a concentração em função do tempo, de t=0 a t=0.5s. Considere h = 0.001;
Observar que os cálculos serão:
Assim, vamos escrever a equação na forma g(x) = 0, onde x = ni+1.
Solução, por exemplo, usando método de Newton
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Ver exemplo Matlab – Euler implicito e Dica_EDO.MDL
Resgatando método de Newton:
Assim:
E:
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.2.5 Método de Range-Kutta
Os métodos de Range-Kutta são métodos explícitos de passo simples cuja solução é calculada na forma: 
 xi+1 = xi + h e yi+1 = yi + Inclinação.h
Para os métodos de Range-Kutta a inclinação é obtida através das inclinações em vários pontos do intervalo.
Os diferentes tipos de métodos de Range-Kutta são classificados de acordo com a sua ordem e esta ordem identifica o número de pontos utilizados para a determinação da inclinação:
 2ª ordem: usa a inclinação em 2 pontos
 3ª ordem: usa a inclinação em 3 pontos
 4ª ordem: usa a inclinação em 4 pontos
...
 N-ésima ordem: usa a inclinação em N pontos
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.2.5a Método de Range-Kutta de 2ª ordem
A fórmula geral é: yi+1 = yi + (c1K1 + c2K2).h Eq - 1
Com: K1 = f(xi,yi) e K2 = f(xi+a2h, yi+b21K1h)
onde c1, c2, a2, b21 são constantes
OBS: estas constantes variam de acordo com o método de 2ª ordem específico
Por exemplo, o método de Euler modificado e o método do Ponto Central são duas versões do método de Range-Kutta de 2ª ordem.
Método de Euler modificado na forma de uma Range-Kutta de 2ª ordem
Use c1 = c2 = ½ e a2 = b21 = 1
Substituindo na Eq-1: com
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Método do Ponto Central na forma de uma Range-Kutta de 2ª ordem
A fórmula geral é: yi+1 = yi + (c1K1 + c2K2).h Eq - 1
Com: K1 = f(xi,yi) e K2 = f(xi+a2h, yi+b21K1h)
Use c1 = 0; c2 = 1 e a2 = b21 = ½ 
Substituindo na Eq-1: com
Algoritmo:
Defina c1, c2, a2, b21 
Calcule K1= f(xi,yi) 
Calcule K2 = f(xi+a2h, yi+b21K1h)
Calcule a inclinação final 
Calcule a solução numérica:
Ver exemplo Matlab – rk2.m
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Dedução do método de Range-Kutta de 2ª ordem
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.2.5b Método de Range-Kutta de 3ª ordem
A fórmula geral é: yi+1 = yi + (c1K1 + c2K2 +c3K3).h 
Com: K1 = f(xi,yi) ; K2 = f(xi+a2h, yi+b21K1h) e K3 = f(xi+a3h,yi + b31K1h +b32K2h)
onde c1, c2, c3, a2, a3, b21, b31, b32 são constantes
OBS: estas constantes variam de acordo com o método de 3ª ordem específico
Para o método de Range-Kutta de 3ª ordem clássico:
c1 = c3 = 1/6 ; c2 = 4/6; a2 = b21 = 1/2 ; a3 = 1; b31 = -1; b32 = 2
Com estas constantes, o método de RK de 3ª ordem clássico fica:
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
A tabela ilustra constantes para alguns métodos de Range-Kutta de 3ª ordem.
Algoritmo:
Defina c1, c2, c3, a2, a3, b21 , b31, b32 
Calcule K1= f(xi,yi) 
Calcule K2 = f(xi+a2h, yi+b21K1h)
Calcule K3 = f(xi+a3h, yi+b31K1h+b32K2h)
Calcule a inclinação final 
Calcule a solução numérica:
Ver exemplo Matlab – rk3.m
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.2.5c Método de Range-Kutta de 4ª ordem
A fórmula geral é: yi+1 = yi + (c1K1 + c2K2 +c3K3 + c4K4).h 
Com: 	K1 = f(xi,yi) ; 
	K2 = f(xi+a2h, yi+b21K1h) ; 
	K3 = f(xi+a3h,yi + b31K1h +b32K2h) ;
	K4 = f(xi+a4h, yi + b41K1h + b42K2h + b43K3h)
onde c1, c2, c3, c4, a2, a3, a4, b21, b31, b32, b41, b42, b43 são constantes
OBS: estas constantes variam de acordo com o método de 4ª ordem específico
Para o método de Range-Kutta de 4ª ordem clássico:
c1 = c4 = 1/6 ; c2 = c3 = 2/6; a2 = a3 = b21 = b32 = 1/2 ; a4 = b43 = 1; b31=b41=b42 = 0
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Substituindo as constantes, o método de RK de 4ª ordem clássico fica:
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Algoritmo:
Defina c1, c2, c3, c4, a2, a3, a4, b21 , b31, b32, b41, b42, b43 
Calcule K1= f(xi,yi) 
Calcule K2 = f(xi+a2h, yi+b21K1h)
Calcule K3 = f(xi+a3h, yi+b31K1h+b32K2h)
Calcule K4 = f(xi+a4h, yi+b41K1h+b42K2h+b43K3h)
Calcule a inclinação final 
Calcule a solução numérica:
Ver exemplo Matlab – rk4.m
Observe que o erro é pequeno mesmo quando o passo de integração foi grande (h=0,5)
Tal fato, por sua vez, não associa menor esforço computacional já que para cada passo de integração foram realizados muito mais cálculos
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.2.5d Método de Range-Kutta de ordem superior (5ª ordem)
OBS: utilizada apenas quando resultados mais acurados forem necessários
A fórmula geral é:
Com:
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.3 Métodos de passo múltiplo
Nos métodos de passo simples, buscávamos a solução no ponto (xi+1, yi+1) utilizando informações do ponto imediatamente anterior (xi, yi).
Nos métodos de passo múltiplo são utilizadas informações de mais pontos anteriores.
Por exemplo:
Num método de passo múltiplo de 2ª ordem, o ponto (xi+1, yi+1) será calculado com base nos pontos (xi, yi) e (xi-1, yi-1)
Num método de passo múltiplo de 3ª ordem, o ponto (xi+1, yi+1) será calculado com base nos pontos (xi, yi) , (xi-1, yi-1) e (xi-2, yi-2)
ENG D04 – Métodos matemáticos e computacionais na engenhariaLeizer Schnitman Resolução Numérica de EDOs
pg.*
3.3.1 Método de Adams-Bashforth (multipasso explícito)
As fórmulas deste método são deduzidas a partir da integração da equação diferencial ao longo de um intervalo arbitrário [xi, xi+1]
Assim, considere a equação diferencial dada na forma
a integração de xi a xi+1 resulta em:
e a integração é realizada representando-se f(x,y) de forma aproximada por um polinômio que interpola esta função em (xi, yi) e em pontos anteriores.
Ou seja: calculamos yi+1 = yi + inclinação.h . Neste caso, porém, a inclinação é também função de pontos anteriores
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Adams-Bashforth de 2ª ordem
Considere que a função f(x,y) é aproximada por um polinômio de primeira ordem (uma reta). Assim:
Fazendo h = xi – xi-1 e substituindo em obtemos:
Que resulta em:
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Observar que foi dada apenas a condição inicial e o método exige que pelo menos 2 pontos sejam conhecidos.
Solução: calcular o primeiro ponto com método de passo simples
(xi,yi)
(xi+1,yi+1)
(xi-1,yi-1)
Ver exemplo adams_b2.m
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Adams-Bashforth de 3ª ordem
Considere que a função f(x,y) é aproximada por um polinômio de segunda ordem. Assim, o cálculo do próximo ponto (i+1) se baseia no ponto atual e em 2 pontos anteriores
Adams-Bashforth de 4ª ordem
Considere que a função f(x,y) é aproximada por um polinômio de teceira ordem. Assim, o cálculo do próximo ponto (i+1) se baseia no ponto atual e em 3 pontos anteriores
Observar problemas de convergência em função do tamanho do passo
Ver exemplo adams_b3.m
Ver exemplo adams_b4.m
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.3.2 Método de Adams-Moulton (multipasso implícito)
O método se assemelha ao método de Adams-Bashforth mas, como método implícito, além do valor atual e de valores anteriores, o cálculo do ponto (i+1) consideta também a inlinação prevista para este ponto (variável aperece dos 2 lados da equação). 
Adams-Moulton de 2ª ordem
Adams-Moulton de 3ª ordem
Adams-Moulton de 4ª ordem
Dá para resolver?
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.3.3 Método de Predição-correção
Vimos que as fórmulas implícitas usualmente requerem maior esforço para sua resolução. O método de predição correção é uma combinação entre uma fórmula explícita e uma fórmula implícita, favorecendo sua implementação.
Genericamente, calcula-se yi+1 através de uma fórmula explícita de passo simples ou multipasso. Este é chamado de valor predito, ou predição.
O valor predito será então corrrigido, ou seja, utiliza-se uma fórmula implícita para calcular o novo yi+1 , onde o valor predito é utilizado no lado direito da equação.
Esta abordagem permite o aproveitamento de vantagens da fórmula implícita e, ao mesmo tempo, evita as dificuldades associadas a solução direta da equação implícita.
Observe que pode-se aplicar o corretor tantas vezes quanto desejadas, de modo a se obter valores de yi+1 tão precisos quanto se queira.
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Para exemplificar o método de predição-correção, vamos utilizar o método de Euler explícito com preditor e o de Euler modificado como corretor
Euler_implicito.m (usou método de Newton)
Algoritmo para o exemplo proposto:
Calcule yi+1 usando o método de Euler_explícito
Inicio do Loop
Calcule o valor corrigido de yi+1 usando o valor predito e o método de Euler_implícito
Avalie a variação entre o valor predito e o valor corrigido
Se a variação for pequena, saia do loop
Caso a variaçào ainda seja grande, faça uma nova “correção”
Fim do Loop
Euler_explicito.m
Ver exemplo Predicao_Correcao.m
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
3.4 Recursos Matlab
OBS: a forma de comando é a mesma para todas as funções
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Comando: [x,y] = NomeDaFuncao (‘EqDif’, X,Yini)
EqDif é um arquivo .m que contém a equação diferencial a ser resolvida
X representa o domínio da variável independente. Por padrão, deve ter dimensão 2 para representar o intervalo de interesse. Pode, porém, ter dimensão maior do que 2, representando pontos intermediários. Neste caso, o cálculo será feito apenas nestes pontos.
Yini representa o valor inicial da função
[x,y] são vetores coluna. O espaçamento depende do domínio definido para X
Para os exemplos seguintes vamos assumir a função EqDif:
function dydt=EqDif(x,y);
dydx = -1.2*y + 7*exp(-0.3*x);
Ver EdoMatlab.m
OBS: existem mais recursos/parâmetros oferecidos pela Matlab e que podem ser explorados
ENG D04 – Métodos matemáticos e computacionais na engenharia Leizer Schnitman Resolução Numérica de EDOs
pg.*
Avaliar:
Sistemas de equações diferenciais
Solução de sistemas de equações diferenciais não-lineares
Conceito de “rigidez” (numa única EDO ou num sistema)
Métodos de RK adaptativos (ajuste do tamanho do passo)
Planos de fase (visualiza estabilidade da solução)
Problemas de valores de contorno (PVC)