Baixe o app para aproveitar ainda mais
Prévia do material em texto
Universidade Federal do Paraná � UFPR Programa de Pós-Graduação em Engenharia Química � PPGEQ Métodos Numéricos em Engenharia Química Prof. Éliton Fontana 2018/1 Conteúdo 1. Introdução 3 1.1. Classi�cação das Equações Diferenciais . . . . . . . . . . . . . . . . . . . . . 4 1.1.1. Equações Diferenciais Ordinárias (EDO) e Parciais (EDP) . . . . . . 4 1.1.2. Ordem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.1.3. Linearidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.4. Homogeneidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.1.5. Coe�cientes constantes e variáveis . . . . . . . . . . . . . . . . . . . . 6 I. Métodos Numéricos para Análise de Equações Diferenciais 7 2. Métodos Numéricos Para PVI's 8 2.1. Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.1. Erro Associado ao Método de Euler . . . . . . . . . . . . . . . . . . . 11 2.1.2. Método de Euler para Sistemas de EDO's . . . . . . . . . . . . . . . 13 2.2. Derivações do Método de Euler . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.2.1. Método de Euler Aprimorado (Fórmula de Heun) . . . . . . . . . . . 16 2.2.2. Método de Euler Implícito (Inverso) . . . . . . . . . . . . . . . . . . . 18 2.2.3. Regra do Ponto Médio . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3. Métodos de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.3.1. Método de Runge-Kutta de 2a Ordem . . . . . . . . . . . . . . . . . 22 2.3.2. Método de Runge-Kutta de 4a Ordem . . . . . . . . . . . . . . . . . 24 2.3.3. Método de Runge-Kutta para Sistemas . . . . . . . . . . . . . . . . . 25 2.3.4. Métodos de Runge-Kutta Adaptativos . . . . . . . . . . . . . . . . . 28 2.4. Métodos de Passos Múltiplos . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4.1. Método de Adams-Bashforth . . . . . . . . . . . . . . . . . . . . . . . 33 2.4.2. Método de Adams-Moulton . . . . . . . . . . . . . . . . . . . . . . . 34 2 2.5. Estabilidade dos Métodos de Resolução de EDO's . . . . . . . . . . . . . . . 37 2.6. Análise de Estabilidade de Métodos de Passo Único . . . . . . . . . . . . . . 38 2.6.1. Equação Modelo: Decaimento de Primeira Ordem . . . . . . . . . . . 39 2.7. Problemas Rígidos (Sti�) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 2.7.1. Sistemas de Equações Diferenciais Rígidas Lineares . . . . . . . . . . 43 2.7.2. Sistemas Não-Lineares . . . . . . . . . . . . . . . . . . . . . . . . . . 44 3. Métodos Numéricos para Problemas de Valor de Contorno 46 3.1. Estratégias de Solução de PVC's . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2. Aproximações por Diferenças Finitas . . . . . . . . . . . . . . . . . . . . . . 48 3.2.1. Aproximação da Derivada Primeira . . . . . . . . . . . . . . . . . . . 49 3.2.2. Aproximação da Derivada Segunda . . . . . . . . . . . . . . . . . . . 52 3.2.3. Discretização das Condições de Contorno . . . . . . . . . . . . . . . . 53 3.2.4. Discretização de PVC's utilizando Diferenças Finitas . . . . . . . . . 54 4. Método das Linhas 60 3 1. Introdução A maioria dos sistemas de interesse na engenharia envolvem mais de uma variável de- pendente que são função da mesma variável independente, como por exemplo a variação na concentração de diversas espécies químicas em um reator ao longo do tempo ou a variação nas três componentes do vetor velocidade ao longo de uma dada direção espacial. De forma geral, não é possível obter a solução para uma única variável independentemente, pois usu- almente existe uma dependência de uma sobre a outra. Isto faz com que os modelos gerem sistemas de equações que devem ser resolvidas ao mesmo tempo, podendo estas equações serem algébricas, diferencias, integrais ou uma mistura delas. Além disso, uma equação diferencia de ordem n pode sempre ser transforma em um sis- tema de n equações diferenciais de primeira ordem. Esta abordagem é particularmente útil para obter soluções numéricas de Problemas de Valor Inicial (PVI's), já que os métodos de resolução (Euler, Runge-Kutta, etc.) são baseados na aproximação da derivada primeira e portanto só podem ser utilizados para resolver equações de primeira ordem. No momento, o objetivo será a análise de sistemas de equações diferenciais ordinárias (EDO's), ou seja, equações que possuem somente uma variável independente. Esta análise pode ser conduzida de forma quantitativa ou qualitativa, sendo que cada forma possui suas vantagens e desvantagens. A análise quantitativa se refere à resolução direta das equações, seja por via analítica ou numérica, de forma a se obter como as variáveis dependentes va- riam em função da independente. Esta forma de análise é útil quando se deseja conhecer esta dependência para uma condição especí�ca, porém trás pouca informação sobre o comporta- mento global do sistema e como uma variável in�uencia as demais. A análise qualitativa, por sua vez, consiste em examinar o comportamento geral de todas as famílias de solução do sistema, o que permite obter informações valiosas sobre a estabilidade das soluções. O material a seguir será dividido em duas partes. Na primeira parte, serão apresentados métodos numéricos que podem ser aplicados na análise de equações diferenciais ordinárias, além do Método das Linhas, que pode ser utilizado para transformar equações diferenciais 4 parciais em um sistema de EDO's. Na segunda parte serão apresentados métodos qualitativos de análise e uma introdução à análise de estabilidade de equações não-lineares. Ao longo do texto serão feitas referências às principais características das equações diferenciais, pois a forma de análise pode depender de determinadas classi�cações. Para facilitar, a seguir será apresentada uma breve revisão da classi�cação das equações diferenciais . 1.1. Classi�cação das Equações Diferenciais As equações diferenciais são classi�cadas de acordo com suas características. Esta classi- �cação é essencial para a determinação dos métodos de análise mais adequados para cada equação. 1.1.1. Equações Diferenciais Ordinárias (EDO) e Parciais (EDP) Quando a função desconhecida depende somente de uma variável independente a equação é chamada de Equação Diferencial Ordinária (EDO). Por exemplo: 4 d2y dt2 + sin(t) dy dt = et De forma genérica, uma EDO pode ser expressa como: f(t, y, y′, y′′, . . . , y(n)) = 0 Caso a variável dependente dependa de duas ou mais variáveis, a equação é classi�cada como Equação Diferencial Parcial (EDP). Por exemplo: ∂T ∂t = α2 ∂2T ∂x2 1.1.2. Ordem A ordem de uma equação diferencial é a ordem da derivada de maior ordem que aparece na equação. Por exemplo: d2y dt2 = 0 segunda ordem dy dt + y3 = 0 primeira ordem A ordem da equação não está relacionada com a maior potência na qual a variável está elevada. 5 1.1.3. Linearidade Uma equação algébrica é dita linear quando pode ser escrita da forma: a1x1 + a2x2 + a3x3 + . . . anxn + b = 0 ou seja, quando as variáveis x1, x2, . . . xn não aparecem em termos não-lineares. De forma semelhante, a equação diferencial ordinária F (t, y, y′, y′′, . . . , y(n)) = 0 é dita linear se F é uma função linear em relação à variável dependente e suas derivadas (y, y′, y′′, . . . , y(n)). Isto implica que uma EDO linear pode ser escrita da forma: an(t) dny dtn + an−1(t) d(n−1)y dt(n−1) + . . .+ a0(t)y = g(t) ou seja, a variável y ou suas derivadas não estão presentes em termos não-lineares, como por exemplo produto entre a variável e as derivadas, yn com n 6= 1 e funções não lineares contendo a variável y (funções trigonométricas, exponencial, etc.). A variável independente (t) pode estar presente em termos não-lineares. Exemplos de equações lineares: d2y dt2 + t dy dt + y = 0 dy dt + t2y = sin(t) dy dt − yet = 0 Quando as equações não podem ser expressas da forma apresentada anteriormente são consideradas não-lineares. As equações diferenciaisnão-lineares possuem em geral uma re- solução muito mais complexa e por isso são mais difíceis de serem avaliadas analiticamente. Exemplos de equações não-lineares: dy dt + y2 = 0 y dy dt = t d2y dt2 = sin(y) dy dt + tey = 0 De forma similar, uma EDP é dita não-linear quando possui algum termo não linear envolvendo qualquer variável dependente. As EDO's lineares ainda podem ser divididas em outras categorias que são utilizadas para facilitar a escolha dos métodos de solução: 1.1.4. Homogeneidade Uma equação diferencial do tipo F (t, y, y′, y′′, . . . , y(n)) = 0 6 é homogênea quando o termo associado somente à variável independente é nulo. Na forma apresentada anteriormente para as EDO's lineares, a equação é homogênea quando g(t) = 0. Caso algum termo diferente de zero não estiver multiplicado pela variável dependente (y) ou alguma de suas derivadas, a equação é dita não-homogênea; 1.1.5. Coe�cientes constantes e variáveis Esta classi�cação costuma ser empregada somente para EDO's. Na expressão para uma EDO genérica vista anteriormente, os coe�cientes an(t), an−1(t)...a0(t) podem ser funções conhecidas da variável independente t. No entanto, caso estes coe�cientes sejam constantes (não dependam de t) a equação é conhecida como EDO com coe�cientes constantes, caso contrário é chamada de EDO com coe�cientes variáveis. Uma classe de equações muito importante nas ciências exatas são as equações autônomas. Estas equações surgem quando a variável independente não aparece de forma explícita na equação. 7 Parte I. Métodos Numéricos para Análise de Equações Diferenciais 8 2. Métodos Numéricos Para PVI's 2.1. Método de Euler A maior parte das equações diferenciais de interesse na engenharia não possuem solução analítica conhecida. Além disso, em muitos casos as soluções analíticas co-nhecidas são difíceis de serem utilizadas, por exemplo quando envolvem séries in�nitas ou integrais sem resolução analítica. Nestes casos é mais conveniente a utilização de ferramentas numéricas para a resolução das equações. Existem basicamente três formas de analisar equações diferenciais: através de métodos analíticos, métodos qualitativos e métodos numéricos. Os métodos analíticos permite esta- belecer uma relação direta entre as variáveis dependentes e independentes através de funções válidas em um determinado intervalo. Os métodos qualitativos envolvem a análise do com- portamento geral da equação sem necessariamente buscar uma solução, como por exemplo através da construção do campo de direções. Os métodos de resolução numérica permitem que a solução seja estimada para uma dada condição inicial e um conjunto de parâmetros especí�cos. A ideia geral é de que é possível obter aproximações da solução de uma determinada EDO avançando em pontos que estão a uma distância ∆t da condição inicial. Diferentemente dos métodos analíticos, a solução obtida através de métodos numéricos é válida somente para um conjunto de parâmetros especí�cos e para uma dada condição inicial. Caso esta condição inicial ou algum parâmetro sejam alterados, deve-se resolver o problema novamente. A determinação de qual forma de análise deve ser utilizada irá depender tanto das caracte- rísticas das equações diferenciais que estão sendo avaliadas quanto de que tipo de informação está se buscando. Assim como para o caso da resolução de sistemas de equações algébricas, existem diversos métodos que podem ser utilizados para a resolução de sistemas de equações diferenciais. Inicialmente será apresentado o método de Euler. Este método costuma ser pouco e�ciente 9 para a resolução de sistemas complexos, porém diversos conceitos fundamentais dos métodos de resolução podem ser mais facilmente apresentados através deste método. O método de Euler é o método mais simples para a aproximação da solução de EDO's. Considere o seguinte problema de valor inicial de primeira ordem: y′ = dy dt = f(t, y) y(t0) = y0 A solução deste PVI passa obrigatoriamente pelo ponto (t0, y0). O objetivo do método de Euler é estimar o valor da solução y(t) em um ponto com uma distância ∆t do ponto inicial t0. Este ponto passa a ser chamado de t1 = t0 + ∆t, de modo que o do método de Euler é determina o valor da solução em t1, y1 = y(t1). A partir deste valor, pode-se estimar o valor da função em um ponto t2 distante ∆t do ponto t1 e continuar com este procedimento até se obter uma aproximação para a solução por quantos pontos forem necessários. Esta classe de métodos são chamados de métodos passo-a-passo. Para avaliar o termo y1 = y(t0 + ∆t), pode-se aplicar uma expansão em séries de Taylor em torno do ponto t0: y(t1) = y(t0 + ∆t) = y0 + ∆ty ′(t0) + ∆t2 2! y′′(t0) + ∆t3 3! y′′′(t0) + . . . Considerando que o espaçamento ∆t seja muito pequeno, os termos ∆t2, ∆t3, etc. podem ser desprezados. Com isso: y(t0 + ∆t) = y1 ≈ y0 + ∆ty′(t0) Considerando também que y′ = f(t, y), a relação anterior pode ser reescrita como: y1 = y0 + ∆tf(t0, y0) De forma geral, a relação anterior pode ser utilizada para estimar a solução para qualquer valor yn+1 com base no valor para yn: yn+1 = yn + ∆tf(tn, yn) Como o valor no novo ponto yn+1 é calculado com base somente em valores no ponto anterior (já conhecidos), este método é dito explícito. Estrutura de um código para resolver uma EDO dy/dt = f(t, y) utilizando o método de Euler: Passo 1: De�nir f(t, y); Passo 2: De�nir os valores iniciais t[1] = t0, y[1] = y0; 10 Passo 3: De�nir o tamanho do passo ∆t e o número de passos n; Passo 4: Para i = 1 até i = n calcular: k1 = f(t[i], y[i]) y[i+ 1] = y[i] + ∆t× k1 t[i+ 1] = t[i] + ∆t (2.1) Passo 5: Plotar os resultados. Exemplo 01): Utilizando o método de Euler, obtenha a solução aproximada até t = 0.5 para o seguinte PVI, utilizando ∆t = 0.1: dy dt = −2t− y y(0) = −1 Neste caso f(t, y) = −2t− y. No ponto inicial temos: t0 = 0 y0 = −1 Como o passo ∆t é constante, os valores de tn são automaticamente de�nidos, t1 = t0 + ∆t = 0.1, t2 = t1 + ∆t = 0.2, . . . Avaliando y1 = y(t1) = y(0.1): y1 = y0 + ∆tf(t0, y0) = −1 + 0.1(−2× 0− (−1)) = −0.9 Repetindo para os próximos passos: y2 = y1 + ∆tf(t1, y1) = −0.9 + 0.1(−2× 0.1− (−0.9)) = −0.83 y3 = y2 + ∆tf(t2, y2) = −0.83 + 0.1(−2× 0.2− (−0.83)) = −0.787 y4 = y3 + ∆tf(t3, y3) = −0.787 + 0.1(−2× 0.3− (−0.787)) = −0.7683 y5 = y4 + ∆tf(t4, y4) = −0.7683 + 0.1(−2× 0.4− (−0.7683)) = −0.77147 Para comparação, a solução exata da equação diferencial é: y(t) = −3e−t − 2t+ 2 sendo que para os pontos avaliados a solução exata é: y0 = −1 y1 = −0.9145 y2 = −0.85619 y3 = −0.8224 y4 = −0.8109 y5 = −0.81959 11 Exemplo 02): Utilizando o método de Euler, obtenha a solução aproximada para t = 1 para o seguinte PVI, utilizando ∆t = 0.2: dy dt = 2t2y2 y(0) = 1 Com base nos dados fornecidos, t0 = 0 e y0 = 1. Como ∆t = 0.2, temos que t1 = 0.2, t2 = 0.4, t3 = 0.6, t4 = 0.8 e t5 = 1. Avaliando os termos intermediários: y1 = y0 + 0.2 · (2 · t20y20) = 1 + 0.2 · (2 · 02 · 12) = 1 y2 = y1 + 0.2 · (2 · t21y21) = 1 + 0.2 · (2 · 0.22 · 12) = 1.016 y3 = y2 + 0.2 · (2 · t22y22) = 1 + 0.2 · (2 · 0.42 · 1.0162) = 1.0821 y4 = y3 + 0.2 · (2 · t23y23) = 1 + 0.2 · (2 · 0.62 · 1.08212) = 1.2507 y5 = y4 + 0.2 · (2 · t24y24) = 1 + 0.2 · (2 · 0.82 · 1.25072) = 1.6511 A solução exata em t = 1 é y(1) = 3, portanto a solução obtida está com um erro signi�cativo. Por este exemplo, �ca claro que é fundamental avaliar os erros associados à utilização do método de Euler. 2.1.1. Erro Associado ao Método de Euler A utilização de métodos numéricos sempre leva à obtenção de soluções aproximadas. Uma importante propriedade dos métodos numéricos é a convergência. Um método é dito convergente quando a solução obtida tende à solução exata quando o espaçamento ∆t→ 0. Na resolução de EDO's pelo método de Euler (ou por qualquer outro método), existem duas fontes de erros: - Erros de truncamento: erros associados com o truncamento da expansão em série de Taylor no segundo termo. Como a determinaçãodo ponto yn+1 depende do valor de yn, o erro de truncamento pode aumentar rapidamente ao longo da resolução. Para reduzir os erros de truncamento, pode-se reduzir o passo de tempo ∆t. - Erros de arredondamento: erro associado à precisão dos valores numéricos utilizados (número de dígitos signi�cativos: 6-9 para precisão única e 15-17 para precisão dupla). Quanto maior for o número de operações necessárias, maior será o erro de arredondamento. Assim, uma maneira de reduzir este erro é aumentar o passo de tempo, de modo que menos passos precisam ser resolvidos para atingir o tempo desejado. 12 Dessa forma, conforme o passo de tempo aumenta, o erro de arredondamento diminui e o de truncamento aumenta. Por isso, existe um ponto ótimo onde o erro total é mínimo, como ilustrado na �gura a seguir. Usualmente os erros de arredondamento são menos signi�cativos, por isso a tendência é que passos pequenos gerem melhores resultados. Como visto anteriormente, pode-se obter o valor de uma função y(t) em um ponto tn+1 com base no valor em tn através de uma expansão em séries de Taylor: yn+1 = yn + ∆t y ′(tn) + ∆t2 2! y′′(tn) + ∆t2 3! y′′′(tn) + . . . A diferença entre o valor exato de yn+1 e o valor aproximado obtido com o uso do método de Euler corresponde ao erro local de truncamento (para um determinado valor de n): en+1 = ∆t2 2! y′′(tn) + ∆t2 3! y′′′(tn) + . . . Nesta forma, o erro ainda representa uma série in�nita. Para evitar este problema, pode-se utilizar a fórmula de Taylor com resto de Lagrange: yn+1 = yn + ∆t y ′(tn) + ∆t2 2! y′′(c) onde c ∈ [tn, tn+1]. Este resultado é uma extensão do teorema do valor médio, que diz que dada uma função contínua f de�nida num intervalo fechado [a,b] e diferenciável em (a,b), existe algum ponto c em (a,b) tal que : f ′(c) = f(b)− f(a) b− a Utilizando a relação anterior, o erro en+1 pode ser avaliado como: en+1 = y′′(c) 2 ∆t2 13 Como y′′(c) é uma constante (valor �nito), temos que: en+1 = O(∆t 2) ou seja, o erro local é da ordem de ∆t2. Como cada passo gera um erro local, quanto mais passos forem utilizados, maior será o acúmulo de erros. Por isso, quando o valor yn+1 estiver sendo determinado, o erro associado à determinação do valor de yn também deve ser considerado. Para n passos, o erro global En pode ser avaliado como: En+1 = en+1 · n = O(h2) · n = O(h2) n · h h = O(h)tn onde h = ∆t. Novamente, como tn é um escalar: En+1 = O(h) Assim, o erro global no método de Euler é da ordem do tamanho do passo utilizado. Esta relação mostra que o erro tende a zero conforme h→ 0 e portanto o método é convergente. De forma geral, um dado método é convergente quando: En+1 = O(h p) p > 0 Neste caso, o método possui ordem p, portanto, o método de Euler é um método de primeira ordem. Métodos de ordem superior convergem mais rapidamente para a solução exata (não exigem passos tão pequenos), sendo por isso mais indicados que os métodos de primeira ordem. A utilização de um passo muito pequeno (tendendo a zero) é impraticável, pois exigiria um tempo computacional demasiadamente alto, além de fazer com que os erros de arredon- damento aumentassem rapidamente. Um procedimento simples que pode na maioria dos casos ser utilizado para determinar se o passo utilizado é adequado é resolver o problema com valores de ∆t gradativamente menores. A partir de um determinado ponto as soluções obtidas serão muito parecidas, sendo que uma maior redução em ∆t a partir deste ponto não irá reduzir o erro global de forma signi�cativa e irá aumentar o erro de arredondamento. 2.1.2. Método de Euler para Sistemas de EDO's O método de Euler explícito pode ser estendido para aplicação em sistemas de EDO's de primeira ordem. Considere o seguinte sistema de PVI's: dx dt = f(t, x, y) x(0) = x0 14 dy dt = g(t, x, y) y(0) = y0 De forma semelhante ao realizado para uma única equação, pode-se partir das condições iniciais e ir avançando ao longo de t com base nos valores já conhecidos: xn+1 = xn + ∆t · f(tn, xn, yn) yn+1 = yn + ∆t · g(tn, xn, yn) Generalizando para um sistema de m equações da forma: dy1 dt = f1(t, y 1, y2, y3, . . . , ym) dy2 dt = f2(t, y 1, y2, y3, . . . , ym) dy3 dt = f3(t, y 1, y2, y3, . . . , ym) ... dym dt = fm(t, y 1, y2, y3, . . . , ym) y1(t0) = y 1 0 y2(t0) = y 2 0 y3(t0) = y 3 0 ... ym(t0) = y m 0 (2.2) A utilização do método de Euler explícito implica em: y1n+1 y2n+1 y3n+1 ... ymn+1 = y1n y2n y3n ... ymn + ∆t f1(tn, y 1 n, y 2 n, . . . , y m n ) f2(tn, y 1 n, y 2 n, . . . , y m n ) f3(tn, y 1 n, y 2 n, . . . , y m n ) ... fm(tn, y 1 n, y 2 n, . . . , y m n ) (2.3) A aplicação do método para resolução de sistemas também permite que ele seja utilizado para a resolução de equações de ordem maior que 1. Por exemplo, considere o seguinte PVI: d2y dx2 = f ( x, y, dy dx ) y(0) = y0 dy dx ∣∣∣ x=0 = y1 Esta equação pode ser transformada em um sistema de primeira ordem através da de�nição de uma variável adicional u: u = dy dx → d 2y dx2 = du dx Com isso, o PVI pode ser reescrito como: dy dx = u y(0) = y0 15 du dx = f (x, y, u) u(0) = y1 Este procedimento pode ser aplicado para escrever uma EDO de qualquer ordem como um sistema de EDO's de primeira ordem. Exemplo 03: Utilizando o método de Euler, obtenha a solução aproximada para os três primeiros passos de tempo para o seguinte sistema, considerando ∆t = 0.1: dx dt = y x(0) = 0 dy dt = x √ y + 3y y(0) = 2 Desse modo, f(t, x, y) = y e g(t, x, y) = x √ y + 3y. Neste caso, temos que t0 = 0, x0 = 0 e y0 = 2. Para avaliar os próximos passos, deve-se determinar todas as variáveis em t1, na sequência todas as variáveis em t2 e assim sucessi- vamente: x1 = x0 + ∆t · f(t0, x0, y0) = 0 + 0.1 · (2) = 0.2 y1 = y0 + ∆t · g(t0, x0, y0) = 2 + 0.1 · (0 · √ 2 + 3 · 2) = 2.6 x2 = x1 + ∆t · f(t1, x1, y1) = 0.2 + 0.1 · (2.6) = 0.46 y2 = y1 + ∆t · g(t1, x1, y1) = 2.6 + 0.1 · (0.2 · √ 2.6 + 3 · 2.6) = 3.4122 x3 = x2 + ∆t · f(t2, x2, y2) = 0.46 + 0.1 · (3.4122) = 0.80122 y3 = y2 + ∆t · g(t2, x2, y2) = 3.4122 + 0.1 · (0.46 · √ 3.4122 + 3 · 3.4122) = 4.520896 16 2.2. Derivações do Método de Euler A utilização do método de Euler para a resolução de problemas típicos da engenharia usualmente implica no uso de passos de tempo muito pequenos para garantir a convergência. Normalmente, melhores resultados são obtidos com o uso de métodos de maior ordem ou que ao menos possuam características de convergência mais atrativas. Como visto na aula passada, o método de Euler possui diversas limitações que restringem o seu uso a problemas muito especí�cos. Dentre as principais causas destas limitações pode- se citar os termos desprezados na expansão em série de Taylor, a necessidade de utilizar passos de tempo muito pequenos e o fato de que o método considera que o valor da derivada y′(ti) é constante em todo o intervalo (ti, ti+1). O método de Euler serviu como base para o desenvolvimento de métodos gradativamente mais complexos e com melhor desempenho. O método de Euler e outros métodos explícitos de maior ordem fazem parte de uma família de métodos numéricos mais abrangentes chamados de métodos de Runge-Kutta (RK). O método de Euler corresponde ao método RK de primeira ordem. Além dos métodos RK, outras abordagens podem ser empregadas, como por exemplo a utilização de métodos implícitos e os métodos de passos múltiplos. A seguir serão apresentadas três modi�cações do método de Euler que pertencem a estas categorias. Assim como para o método de Euler, o objetivo continua sendo buscar uma solução apro- ximada para o PVI: dy dt = f(t, y) y(t0) = y0 2.2.1. Método de Euler Aprimorado (Fórmula de Heun) Este método faz parte da categoria de passos múltiplos, onde os valores da derivada em mais de um ponto são utilizados para determinar o próximo ponto. A fórmula de Heun é uma modi�cação simples do método de Euler e consisteem uma abordagem preditiva-corretiva. O passo inicial consiste determinar um valor preditor com base no método de Euler: y∗n+1 = yn + ∆t · f(tn, yn) A partir deste valor, é calculada uma etapa corretora com base no valor médio da função avaliada em yn e no valor preditor y∗n+1: yn+1 = yn + ∆t 2 · (f(tn, yn) + f(tn+1, y∗n+1)) 17 Esta fórmula representa uma correção do valor estimado anteriormente y∗n+1. Para implementar computacionalmente este método, pode-se utilizar a mesma estrutura geral apresentada para o método de Euler, sendo necessário somente modi�car o passo 4: Passo 4: Para i = 1 até i = n calcular: k1 = f(t[i], y[i]) yp[i+ 1] = y[i] + ∆t× k1 t[i+ 1] = t[i] + ∆t k2 = f(t[i+ 1], yp[i+ 1]) y[i+ 1] = y[i] + (∆t/2)(k1 + k2) Exemplo 01) Repita o exemplo visto anteriormente utilizando a fórmula de Euler apri- morada. Lembrando do PVI (∆t = 0.2): dy dt = 2t2y2 y(0) = 1 A resolução com o método de Euler modi�cado envolve mais etapas: - Passo 1: k1 = f(t0, y0) = 0 yp1 = y0 + ∆t · k1 = 1 t1 = t0 + ∆t = 0.2 k2 = f(t1, yp1) = 0.08 y1 = y0 + (∆t/2) · (k1 + k2) = 1.008 - Passo 2: k1 = f(t1, y1) = 0.08128 yp2 = y1 + ∆t · k1 = 1.0243 t2 = t1 + ∆t = 0.4 k2 = f(t2, yp2) = 0.3357 y2 = y1 + (∆t/2) · (k1 + k2) = 1.0497 - Passo 3: k1 = f(t2, y2) = 0.35259 yp3 = y2 + ∆t · k1 = 1.1202 t3 = t2 + ∆t = 0.6 k2 = f(t3, yp3) = 0.90349 y3 = y2 + (∆t/2) · (k1 + k2) = 1.1753 - Passo 4: k1 = f(t3, y3) = 0.99456 yp4 = y3 + ∆t · k1 = 1.3742 t4 = t3 + ∆t = 0.8 18 k2 = f(t4, yp4) = 2.4172 y4 = y3 + (∆t/2) · (k1 + k2) = 1.5165 - Passo 5: k1 = f(t4, y4) = 2.94371 yp5 = y4 + ∆t · k1 = 2.1052 t5 = t4 + ∆t = 1 k2 = f(t5, yp5) = 8.8637 y5 = y4 + (∆t/2) · (k1 + k2) = 2.6973 Para comparação, a solução exata para t = 0.8 é y = 1.5182 e para t = 1 é y = 3. 2.2.2. Método de Euler Implícito (Inverso) O método de Euler explícito utiliza a inclinação da reta tangente no ponto n para estimar o valor da solução no ponto n+ 1: yn+1 = yn + ∆t · (f(tn, yn)) De forma alternativa, poderia-se utilizara a inclinação da reta tangente no próprio ponto n+ 1. Esta formulação dá origem ao método de Eulet implícito ou inverso ou regressivo: yn+1 = yn + ∆t · (f(tn+1, yn+1)) Como a variável yn+1 aparece nos dois lados da igualdade, pode ser necessário utilizar algum método iterativo para a resolução do problema, caso não seja possível explicitar o valor yn+1. Para EDO's lineares sempre é possível isolar yn+1, no entanto para não-lineares isso usualmente não é possível. Esta formulação representa um método implícito, pois a variável yn+1 depende impli- citamente dela mesmo. Assim como a formulação explícita, este método representa uma aproximação de primeira ordem, porém, como será visto nas próximas aulas, os métodos implícitos possuem uma grande vantagem sobre os explícitos em relação à estabilidade do método. Exemplo 02) Resolva o seguinte PVI utilizando o método de Euler implícito, com ∆t = 0.1 até t = 0.3: dy dt = t2 + ty y(0) = 1 19 2.2.3. Regra do Ponto Médio Lembrando da Fórmula de Taylor com o resto de Lagrange: yn+1 = yn + ∆t y ′(tn) + ∆t2 2! y′′(c) Quando os termos de segunda ordem são desconsiderados, obtém-se o método de Euler explícito. Uma maneira de se conseguir uma aproximação mais precisa (com maior ordem) é considerar um maior número de termos na expansão. Por exemplo, caso o termo de terceira ordem também for considerado, temos que: yn+1 = yn + ∆t y ′(tn) + ∆t2 2! y′′(tn) + ∆t3 3! y′′′(c) onde novamente c ∈ (tn, tn+1). De forma semelhante, pode-se substituir ∆t por −∆t para se obter uma aproximação para yn−1: yn−1 = yn −∆t y′(tn) + ∆t2 2! y′′(tn)− ∆t3 3! y′′′(d) onde d ∈ (tn−1, tn). Fazendo a expansão para yn+1 menos a expansão para yn−1: yn+1 − yn−1 = 2∆ty′(tn) + y′′′(c) + y′′′(d) 3! ∆t3 Novamente, como as derivadas y′′′(c) e y′′′(d) são dois escalares, a relação anterior pode ser expressa como: yn+1 = yn−1 + 2∆ty ′(tn) +O(∆t 3) Desconsiderando os termos de terceira ordem: yn+1 = yn−1 + 2∆ty ′(tn) Esta relação é conhecida como regra do ponto médio e representa uma aproximação de segunda ordem para yn+1. Esta fórmula depende do conhecimento do valor da variável em três pontos, por isso não pode ser utilizara para determinar y1, já que y−1 não existe. Outra relação pode ser utilizada para determinar o ponto y1, e a relação anterior pode ser usada para os demais pontos. Comparando com o método de Euler explícito: yn+1 = yn + ∆ty ′(tn) → y′(tn) = yn+1 − yn ∆t 20 yn+1 = yn−1 + 2∆ty ′(tn) → y′(tn) = yn+1 − yn−1 2∆t Exemplo 07:) Utilize a regra do ponto médio para resolver o seguinte PVI até t = 0.3, utilizando ∆t = 0.1: dy dt = y + 2t+ t2 y(0) = 1 Para avaliar o primeiro ponto, pode-se utilizar alguma formulação distinta. Como a regra do ponto médio é um método de segunda ordem, é recomendável a utilização de outro método de segunda ordem para não haver uma perda muito grande de precisão. Utilizando o método de Euler aprimorado: k1 = f(t0, y0) = 1 yp1 = y0 + ∆t · k1 = 1.1 t1 = t0 + ∆t = 0.1 k2 = f(t1, y1) = 1.31 y1 = y0 + (∆t/2)(k1 + k2) = 1.1155 A partir dos valores conhecidos para y0 e y1, pode-se agora determinar os demais pontos: y2 = y0 + 2∆tf(t1, y1) = 1 + 2 · 0.1 · (1.1155 + 2 · 0.1 + 0.12) = 1.2651 y3 = y1 + 2∆tf(t2, y2) = 1.1155 + 2 · 0.1 · (1.2651 + 2 · 0.2 + 0.22) = 1.4565 21 2.3. Métodos de Runge-Kutta Os métodos de Runge-Kutta (RK) utilizam a aproximação em série de Taylor para obter formulações de alta ordem sem necessitar o cálculo das derivadas de alta ordem. De forma geral, estes métodos podem ser expressos de forma geral como: yi+1 = yi + ∆tφ(ti, yi,∆t) onde a função φ(ti, yi,∆t) é chamada de incremento e pode ser interpretada como uma representação da inclinação no intervalo entre ti e ti+1. Como esta função só depende de pontos já conhecidos, os métodos de Runge-Kutta clássicos são explícitos. Para um método de ordem n, a função incremento costuma ser expressa na seguinte forma: φ(ti, yi,∆t) = a1k1 + a2k2 + a3k3 + . . .+ ankn onde a1, a2, a3, . . . , an são constantes e os termos ki são dados por: k1 = f(ti, yi) k2 = f(ti + p1∆t, yi + q11k1∆t) k3 = f(ti + p2∆t, yi + q21k1∆t+ q22k2∆t) ... kn = f(ti + pn−1∆t, yi + qn−1,1k1∆t+ qn−1,2k2∆t+ . . .+ qn−1,n−1kn−1∆t) onde os termos pi e qij são constantes. Como pode ser observado, a determinação do parâ- metros kn depende os valores ki onde i < n. Por exemplo, para um método de primeira ordem a função incremento é avaliada como: φ(ti, yi,∆t) = a1k1 onde k1 = f(ti, yi). Substituindo na expressão para a forma geral dos métodos de Runge- Kutta: yi+1 = yi + ∆t(a1f(ti, yi)) Os métodos RK estão baseados diretamente na expansão em série de Taylor em torno do ponto yi. Como visto anteriormente, a expansão de primeira ordem resulta em: yi+1 = yi + ∆tf(ti, yi) Comparando os termos semelhantes com a equação anterior, observa-se que a1 = 1. O método de Runge-Kutta de primeira ordem corresponde, de fato, ao método de Euler. 22 2.3.1. Método de Runge-Kutta de 2a Ordem Considere agora que se deseja obter uma aproximação de segunda ordem. Neste caso, temos que: yi+1 = yi + (a1k1 + a2k2)∆t k1 = f(ti, yi) k2 = f(ti + p1∆t, yi + q11k1∆t) Para determinar os valores de a1, a2, p1 e q11 pode-se primeiramente avaliar a expansão em séries de Taylor de segunda ordem em torno de yi é dada por: yi+1 = yi + ∆tf(ti, yi) + df dy ∆t2 2 Considerando que f = f(t, y), a derivada total de f em relação a t é dada por: df dt = ∂f ∂t + ∂f ∂y dy dt Substituindo na expressão anterior: yi+1 = yi + ∆tf(ti, yi) + ( ∂f ∂t + ∂f ∂y dy dt ) ∆t2 2 Utilizando um procedimento semelhante, pode-se utilizar a expansão em séries de Taylor para expandir a própria função f(t, y). Lembrando, para uma função de duas variáveis, a expansão em série de Taylor é dada por: f(x, y) = f(x0, y0) + ∂f ∂x (x− x0) + ∂f ∂y (y − y0)+ 1 2! ( ∂2f ∂x2 (x− x0)2 + ∂2f ∂x∂y (x− x0)(y − y0) + ∂2f ∂y2 (y − y0)2 ) + . . . De forma equivalente, pode-se avaliara função em um ponto qualquer f(x+ a, y+ b) com base no valor de f(x, y). Considerando uma aproximação de primeira ordem: f(x+ a, y + b) = f(x, y) + ∂f ∂x a+ ∂f ∂y b Assim, considerando a de�nição de k2, o termo f(ti+p1∆t, yi+q11k1∆t) pode ser expresso como: k2 = f(ti + p1∆t, yi + q11k1∆t) = f(ti, yi) + p1∆t ∂f ∂t + q11k1∆t ∂f ∂y Substituindo os valores de k1 e k2 na expressão para yi+1 = yi+(a1k1+a2k2)∆t e colocando os termos de mesma ordem em evidência: yi+1 = yi + (a1f(ti, yi) + a2f(ti, yi))∆t+ ( a2p1 ∂f ∂t + a2q11f(ti, yi) ∂f ∂y ) ∆t2 23 Comparando com a equação obtida anteriormente: yi+1 = yi + ∆tf(ti, yi) + ( ∂f ∂t + ∂f ∂y dy dt ) ∆t2 2 Igualando os coe�cientes com os termos de mesma ordem, observa-se que as seguintes relações devem ser satisfeitas: a1 + a1 = 1 a2p1 = 1/2 a2q11 = 1/2 Este sistema possui três equações e quatro variáveis, portanto não existe solução única. No entanto, a partir do momento em que um dos valores é especi�cado, os demais podem ser calculados. Isto implica que existem in�nitas formulações para o método de Runge-Kutta de segunda ordem. Euler Modi�cado (a2 = 0.5) Assumindo que a2 = 0.5, obtém-se que a1 = 0.5 e p1 = q11 = 1. Assim, o método de Runge-Kutta de segunda ordem resulta em: k1 = f(ti, yi) k2 = f(ti + ∆t, yi + k1∆t) yi+1 = yi + (k1 + k2) 2 ∆t Este método corresponde ao método de Euler modi�cado (fórmula de Heun). Regra do Ponto Médio (a2 = 1) Fazendo a2 = 1, temos que a1 = 0 e p1 = q11 = 1/2: k1 = f(ti, yi) k2 = f(ti + ∆t/2, yi + k1∆t/2) yi+1 = yi + k2∆t Esta formulação corresponde ao método de Euler modi�cado com base na regra do ponto médio. 24 Método de Ralston (a2 = 2/3) Fazendo a2 = 2/3, os demais valores são avaliados como a1 = 1/3 e p1 = q11 = 3/4: k1 = f(ti, yi) k2 = f(ti + 3∆t/4, yi + 3k1∆t/4) yi+1 = yi + (k1/3 + 2k2/3)∆t Esta formulação é chamada de Método de Ralston 2.3.2. Método de Runge-Kutta de 4a Ordem Seguindo o mesmo procedimento mostrado anteriormente para uma aproximação de quarta ordem, obtém-se o método RK4. Este é possivelmente o método mais empregado para a resolução de PVI's, devido a sua alta precisão e por ser um método explícito. Neste caso, os parâmetros k1, k2, k3 e k4 são avaliados como: k1 = f(ti, yi) k2 = f(ti + ∆t/2, yi + k1∆t/2) k3 = f(ti + ∆t/2, yi + k2∆t/2) k4 = f(ti+1, yi + k3∆t) Com base nestes valores, o ponto yn+1 é avaliado como: yn+1 = yn + (1/6)(k1 + 2 · k2 + 2 · k3 + k4) A utilização de passos de tempo ∆t constantes para todos os casos possui uma série de desvantagens, especialmente porque em muitos casos a equação pode precisar de valores muito pequenos para manter a precisão em alguns pontos, enquanto que em outros pontos valores maiores poderiam ser utilizados para reduzir o número de cálculos realizados. Muitos softwares utilizam uma abordagem que determina localmente o tamanho necessário para os passos, através da comparação dos resultados obtidos com RK4 e uma formulação de Runge-Kutta de quinta ordem (ex. ode45 do Matlab). Exemplo 08:) Utilize o método RK4 para resolver o seguinte PVI, adotando ∆t = 0.1 até t = 0.3: dy dt = −y + 2t y(0) = 2 25 2.3.3. Método de Runge-Kutta para Sistemas O método de Runge-Kutta pode ser utilizado para a resolução de sistemas de EDO's. Considere, por exemplo, um sistema com duas equações: dy dt = f(y, x, t) y(t0) = y0 dx dt = g(y, x, t) x(t0) = x0 Neste caso os parâmetros k1, . . . , k4 irão depender de qual equação está sendo avaliada: ky1 = f(tn, yn, xn) kx1 = g(tn, yn, xn) ky2 = f(tn + ∆t/2, yn + k y 1∆t/2, xn + k x 1 ∆t/2) kx2 = g(tn + ∆t/2, yn + k y 1∆t/2, xn + k x 1 ∆t/2) ky3 = f(tn + ∆t/2, yn + k y 2∆t/2, xn + k x 2 ∆t/2) kx3 = g(tn + ∆t/2, yn + k y 2∆t/2, xn + k x 2 ∆t/2) ky4 = f(tn + ∆t, yn + k y 3∆t, xn + k x 3 ∆t) kx4 = g(tn + ∆t, yn + k y 3∆t, xn + k x 3 ∆t) Com base nestes valores, a solução no ponto n+ 1 pode ser aproximada: yn+1 = yn + ∆t 6 (ky1 + 2k y 2 + 2k y 3 + k y 4) xn+1 = xn + ∆t 6 (kx1 + 2k x 2 + 2k x 3 + k x 4 ) Considere agora um sistema de m equações: dy1 dt = f1(t, y 1, y2, y3, . . . , ym) dy2 dt = f2(t, y 1, y2, y3, . . . , ym) dy3 dt = f3(t, y 1, y2, y3, . . . , ym) ... dym dt = fm(t, y 1, y2, y3, . . . , ym) y1(t0) = y 1 0 y2(t0) = y 2 0 y3(t0) = y 3 0 ... ym(t0) = y m 0 (2.4) Neste caso, as quantidades k1, k2, etc. passam a ser vetores com m componentes: 26 k1 = f1(xn, y 1 n, y 2 n, . . . , y m n ) f2(xn, y 1 n, y 2 n, . . . , y m n ) ... fm(xn, y 1 n, y 2 n, . . . , y m n ) k2 = ∆t f1(xn + ∆t/2, y 1 n + ∆t · k1[1]/2, y2n + ∆t · k1[2]/2, . . . ymn + ∆t · k1[m]/2) f2(xn + ∆t/2, y 1 n + ∆t · k1[1]/2, y2n + ∆t · k1[2]/2, . . . ymn + ∆t · k1[m]/2) ... fm(xn + ∆t/2, y 1 n + ∆t · k1[1]/2, y2n + ∆t · k1[2]/2, . . . ymn + ∆t · k1[m]/2) k3 = ∆t f1(xn + ∆t/2, y 1 n + ∆t · k2[1]/2, y2n + ∆t · k2[2]/2, . . . ymn + ∆t · k2[m]/2) f2(xn + ∆t/2, y 1 n + ∆t · k2[1]/2, y2n + ∆t · k2[2]/2, . . . ymn + ∆t · k2[m]/2) ... fm(xn + ∆t/2, y 1 n + ∆t · k2[1]/2, y2n + ∆t · k2[2]/2, . . . ymn + ∆t · k2[m]/2) k4 = ∆t f1(xn + ∆t, y 1 n + ∆t · k3[1], y2n + ∆t · k3[2], . . . ymn + ∆t · k3[m]) f2(xn + ∆t, y 1 n + ∆t · k3[1], y2n + ∆t · k3[2], . . . ymn + ∆t · k3[m]) ... fm(xn + ∆t, y 1 n + ∆t · k3[1], y2n + ∆t · k3[2], . . . ymn + ∆t · k3[m]) Com base nestes vetores, a solução para o sistema pode ser avaliada: yin+1 = y i n + ∆t 2 (k1[i] + 2k2[i] + 2k3[i] + k4[i]) Exemplo 01:) Utilize o método RK4 para resolver os dois primeros passos para o seguinte sistema, considerando ∆t = 0.2: dy dt = −7y sin(x+ 2t) y(0) = 1 dx dt = 4x− y2 x(0) = 0 De�nindo: f(t, y, x) = −7y sin(x+ 2t) g(t, y, x) = 4x− y2 t0 = 0 y0 = 1 x0 = 0 Para o primeiro passo de tempo, temos que: t1 = t0 + ∆t = 0.2 27 Avaliando primeiramente k1 para as duas equações: ky1 = f(t0, y0, x0) = (−7 · (1) sin(0 + 2 · 0)) = 0 kx1 = g(t0, y0, x0) = (4 · (0)− 12) = −1 Avaliando agora k2 ky2 = f(t0+∆t/2, y0+∆tk y 1/2, x0+∆tk x 1/2) = (−7·(1) sin(−1·0.2/2+2·0.2/2)) = −0.6988339 kx2 = g(t0 + ∆t/2, y0 + ∆tk y 1/2, x0 + ∆tk x 1/2) = (4 · (−1 · 0.2/2)− (1)2) = −1.4 Para k3: ky3 = ∆tf(t0 + ∆t/2, y0 + ∆tk y 2/2, x0 + ∆tk x 2/2) = (−7 · (1 + 0.2(−0.6988339)/2) sin((0.2 ∗ (−1.4)/2) + 2 ∗ 0.1)) = −0.3904146283 kx3 = ∆tg(t0 + ∆t/2, y0 + ∆tk y 2/2, x0 + ∆tk x 2/2) = (4(0.2 · (−1.4)/2)− (1 + 0.2(−0.6988339)/2)2) = −1.42511690 Para k4: ky4 = f(t1, y0 + ∆tk y 3 , x0 + ∆tk x 3 ) = (−7(1 + 0.2 · (−0.3904146283)) sin((0.2 · (−1.42511690))− 2 · (0.2))) = −0.7403586 kx4 = g(t1, y0 + ∆tk y 3 , x0 + ∆tk x 3 ) = (4(0.2 · (−1.42511690))− (1 + 0.2 · (−0.3904146283))2) = −1.99002461 Agora os pontos y1 e x1 podem ser avaliados: y1 = ∆t 6 (ky1 + 2k y 2 + 2k y 3 + k y 4) = 0.9027048 x1 = ∆t 6 (kx1 + 2k x 2 + 2k x 3 + k x 4 ) = −0.2880086 Repetindo o procedimento para o próximo passo, obtém-se os seguintes coe�cientes: ky1 = −0.706188 kx1 = −1.9669104 k y 2 = −0.6700916 kx2 = −2.6311658 ky3 = −0.285797 kx3 = −2.902888 k y 4 = 0.405631 k x 4 = −4.1892917 Com isso, obtém-se: y2 = 0.828960 x2 = −0.8621522 28 2.3.4. Métodos de Runge-Kutta Adaptativos Como visto anteriormente, o passo de tempo utilizado (∆t) possui uma grande in�uência no erro associado à resolução de equações diferenciais através de métodos numéricos. Além disso, o valor de ∆t irá de�nir quantos passos serão necessários para atingir o tempo �nal desejado. Para reduzir o gasto computacional envolvido, deve-se utilizar o maior passo de tempo possível que garanta que o erro esteja dentro do limite estabelecido. Para um grande número de problemas comuns na engenharia, a resolução pode necessitar de passos muito pequenos em um determinado intervalo de tempo e possibilitar que passos maiores sejam utilizados para os demais intervalos. Por exemplo, quando a solução representa um decaimento exponencial, usualmente nos instantes iniciais ocorre uma rápida variação na solução que tende a estabilizarconforme o tempo avança. Como consequência, passos de tempo pequenos são necessários nos instantes inicias, porém após um determinado tempo os passos podem ser aumentados sem que o erro seja superior ao especi�cado. Para evitar estes inconvenientes, pode-se utilizar um passo de tempo que se adapte a solução. A implementação deste tipo de algoritmo requer que o erro de truncamento local seja de alguma forma estimado ao longo da resolução (para cada passo de tempo). Apesar de isto gerar um gasto computacional extra, normalmente o ganho associado ao uso de passos adaptativos supera em muito este gasto. Existem duas formas duas abordagens distintas para implimentar métodos de Runge- Kutta com passo adaptativo. Uma delas consiste em estimar o erro através da diferença entre duas predições obtidas com passos diferentes usando um método de mesma ordem. A outra abordagem consiste em avaliar o erro local de truncamento através da comparação entre os resultados obtidos com métodos de ordem distintas, utilizado o mesmo passo de tempo. Método de Passo Duplo O método adaptativo mais simples é o método de passo duplo, que consiste em avaliar cada ponto duas vezes, uma vez através de um único passo e outra através de dois passos com a metade do tamanho do primeiro. A diferença entre os dois resultados representa uma estimativa do erro de truncamento local. Por exemplo, quando o ponto yi+1 vai ser calculado com base no ponto yi, primeiramente calcula-se utilizando um passo ∆t e na sequência calcula-se novamente yi+1 através de dois passos com tamanho ∆t/2. 29 Suponha que y∗1 representa o valor calculado em t+∆ com o passo ∆t e y ∗ 2 o valor calculado para o mesmo tempo t+ ∆t mas com o passo ∆t/2, como ilustrado na �gura a seguir. Por exemplo, considere o método de Runge-Kutta de primeira ordem (Euler explícito). Como visto anteriormente, neste caso o erro de truncamento local é da ordem de ∆t2. De�nindo a solução exata em t + ∆t como x, pode-se de�nir a solução exata em termos de y∗1: x = y∗1 + ∆t 2φ+O(∆t3) onde φ é um escalar. De forma semelhante, para y∗2 são necessários dois passos, de modo que: x = y∗2 + 2 ( ∆t 2 )2 φ+O(∆t3) = y∗2 + ∆t2 2 φ+O(∆t3) Desprezando os termos da ordem de ∆t3, igualando as duas expressões temos que: y∗1 + ∆t 2φ = y∗2 + 2 ( ∆t 2 )2 φ De�nindo ∆ = y∗2 − y∗1: ∆ = y∗2 − y∗1 = ∆t2 2 φ Assim, o termo ∆ representa uma estimativa do erro associado com y∗2. Para métodos de maior ordem, o parâmetro ∆ também irá representar uma estimativa do erro associado com y∗2, porém neste caso o erro não será da ordem de ∆t 2 mas irá depender da ordem do método utilizado. Por exemplo, para o método RK4, o parâmetro ∆ será da orde de ∆t5. Se o valor ∆ for inferior ao erro local permitido, pode-se aumentar o passo ∆t, enquanto que se ∆ for maior que o erro permitido, deve-se diminuir ∆t. 30 Considere, por exemplo, que o erro máximo permitido seja emax = 10−6. Se ∆ = |y∗2−y∗1| < 10−6, valida-se o passo e utiliza-se o valor y∗2. Para corrigir o passo de tempo ∆t, diferentes relações podem ser utilizadas, como por exemplo: ∆t = (emax ∆ )α ∆t∗ onde ∆t∗ representa o passo de tempo antigo e ∆t o valor corrigido. Como neste caso emax > ∆, o valor novo obtido será maior que o anterior. O valor de α recomendado para este caso (emax > ∆) é de α = 0.2 Considerando agora que ∆ > emax. Neste caso, o passo é invalidado e deve ser recalculado com um passo menor. A relação utilizada para recalcular o passo é a mesma apresentada anteriormente, porém neste caso o valor corrigido será menor que o anterior e o valor reco- mendado para α é de α = 0.25. Exemplo 02) Utilize o método de Runge-Kutta de quarta ordem juntamente com o método do passo duplo para estimar o valor de y(1.25). Considere que o erro de truncamento local não pode ser superior a 10−3 e como estimativa inicial para o passo de tempo considere ∆t = 1: dy dt = 4e0.8t − 0.5y y(0) = 2 Avaliando o primero passo com ∆t = 1, obtemos y(1) = y∗1 = 6.201. Considerando agora ∆t = 0.5, a resolução de dois passos implica que y(1) = y∗2 = 6.195. Assim: ∆ = y∗2 − y∗1 = 0.006 Como o valor é superior ao erro máximo, deve-se refazer este passo. Recalculando o passo de tempo: ∆t = ( 10−3 0.006 )0.25 1 = 0.638943 Recalculando o primeiro passo de tempo com o novo valor de ∆t, obtém-se que y∗1 = 4.3481 e y∗2 = 4.3475, de modo que ∆ = 6 × 10−4. Como o valor é inferior ao erro máximo, este ponto pode ser validado: t1 = t0 + ∆t = 0.638943 y(t1) = y1 = 4.3475 Para avaliar o próximo passo, primeiramente pode-se avaliar um novo ∆t: ∆t = ( 10−3 6× 10−4 )0.2 · 0.638943 = 0.707672 31 Utilizando este valor, obtém-se qu para t2 = t1 + ∆t = 0.638943 + 0.707672 = 1.346615, y∗1 = 8.4886328 e y ∗ 2 = 8.4869326, de modo que ∆ = 0.0017. Como o valor está acima do erro máximo, deve-se reduzir o passo de tempo: ∆t = ( 10−3 0.0017 )0.25 0.707672 = 0.619755 Recalculado os valores, agora para t2 = t1 +∆t = 0.638943+0.619755 = 1.258698, obtém- se que y∗1 = 7.849350209 e y ∗ 2 = 7.84849122, de modo que ∆ = 8.59 × 10−4. Como o erro está abaixo do especi�cado, pode-se validar o passo. Assim: t2 = 1.258698 y(t2) = y2 = 7.84849122 Como deseja-se saber y(1.25), o valor de t2 já está acima do especifado. Para obter o valor no ponto t = 1.25, pode-se realizar uma interpolação linear entre os valores para t1 e t2: y(1.25) = y1 + (t2 − 1.25) (t2 − t1) (y2 − y1) = 7.79936 Método de Runge-Kutta-Fehlberg Outra alternativa para estimar o erro local de truncamento é avaliar um dado ponto através de métodos com ordem distinta. Por exemplo, pode-se avaliar a solução em um dado ponto y(ti+1) através de um método de segunda ordem (y∗i+1) e de um método de terceira ordem (yi+1). A diferença entre os dois valores (yi+1 − y∗i+1) representa uma estimativa do erro de truncamento neste ponto. Se os dois valores forem muito próximos, isto indica que a redução no passo de tempo não irá diminuir signi�cativamente o erro local de truncamento. A abordagem mais utilizada consiste em avaliar a solução em um ponto através de métodos de Runge-Kutta de quarta e quinta ordem, sendo que as estimativas de quarta ordem (y∗i+1) e de quinta ordem (yi+1) podem ser representadas como: y∗i+1 = yi + ∆t ( 37 378 k1 + 250 621 k3 + 125 594 k4 + 512 1771 k6 ) yi+1 = yi + ∆t ( 2825 27648 k1 + 18575 48384 k3 + 13525 55296 k4 + 277 14336 k5 1 4 + k6 ) onde: k1 = f(ti, yi) 32 k2 = f ( ti + ∆t 5 , yi + ∆t 5 k1 ) k3 = f ( ti + 3∆t 10 , yi + 3∆t 40 k1 + 9∆t 40 k2 ) k4 = f ( ti + 3∆t 5 , yi + 3∆t 10 k1 − 9∆t 10 k2 + 6∆t 5 k3 ) k5 = f ( ti + ∆t, yi − 11∆t 54 k1 + 5∆t 2 k2 − 70∆t 27 k3 + 35∆t 27 k4 ) k6 = f ( ti + 7∆t 8 , yi + 1631∆t 55296 k1 + 175∆t 512 k2 + 575∆t 13824 k3 + 44275∆t 110592 k4 + 253∆t 4096 k5 ) Da mesma forma como para o método anterior, o erro associado com a estimativa do ponto yi+1 é avaliado como: ∆ = yi+1 − y∗i+1 Com base neste valor, utiliza-se o mesmo critério apresentado anteriormente para corrigir o passo de tempo ∆t, ou seja, se o erro for menor que o especi�cado aumenta-se o passo de tempo e se o erro for menor descarta-se o valor obtido e reduz-se o passo de tempo até que o erro esteja dentro do valor tolerável. As relações utilizadas para reduzir ou aumentar o passo são as mesmas apresentadas para o método do passo duplo. 33 2.4. Métodos de Passos Múltiplos Os métodos de passos múltiplos são aqueles onde a variável yn+1 é determinada com base nos valores desta variável em mais de um ponto, como por exemplo no método baseado na regra do ponto médio: yn+1 = yn−1 + 2∆ty ′(tn) → y′(tn) = yn+1 − yn−1 2∆t Os demais métodos, que determinam yn+1 com base somente nos valores de yn, são cha- mados de métodos de passo único. Como visto anteriormente, os métodos de passos múltiplos não iniciam automaticamente, sendo necessário avaliar alguns pontos iniciais através de outro método. Dentre os métodos de passos múltiplos, os mais utilizadossão os da família dos métodos de Adams, como apresentado a seguir. 2.4.1. Método de Adams-Bashforth Considere novamente o seguinte PVI: dy dt = f(t, y) y(t0) = y0 Uma maneira de obter o valor em um ponto yn+1 com base em um valor conhecido yn é integrar a equação desde um ponto tn até um ponto tn+1:∫ tn+1 tn y′(t)dt = ∫ tn+1 tn f(t, y)dt Como o lado esquerdo avalia a integral de uma derivada, podemos escrever a expressão anterior como: y(tn+1)− y(tn) = yn+1 − yn = ∫ tn+1 tn f(t, y)dt Como a função y(t) não é conhecida, não é possível resolver diretamente a integral do lado direito da equação. O método de Adams-Bashforth consiste em substituir a função f(t, y) por um polinômio p(t), permitindo assim a resolução da integral: yn+1 = yn + ∫ tn+1 tn p(t)dt Dependendo da ordem escolhida para o polinômio p(t), diferentes formulações são obtidas. Os coe�cientes associados ao polinômio são determinados com base nos valores já conhecidos 34 para yn, yn−1, yn−2, etc. Caso um polinômio de ordem k for utilizado, é necessário conhecer a solução em k + 1 pontos. Por exemplo, caso um polinômio de primeira ordem for empregado: p(t) = At+B é necessário conhecer o valor da função f(x, y) em dois pontos, ou seja, é necessário deter- minar dois valores (yn e yn−1) para encontar as constantes A e B. Com isso: Atn +B = f(tn, yn) = fn Atn−1 +B = A(tn −∆t) +B = f(tn−1, yn−1) = fn−1 Resolvendo para A e B: A = fn − fn−1 ∆t B = fn−1tn − fntn−1 ∆t Substituindo p(t) = At+B na expressão anterior e resolvendo a integral: yn+1 = yn + ∫ tn+1 tn (At+B)dt = yn + A 2 (t2n+1 − t2n) +B(tn+1 − tn) Substituindo as expressões obtidas para A e B e simpli�cando o resultado obtido: yn+1 = yn + 3 2 ∆t · fn − 1 2 ∆t · fn−1 Esta relação é conhecida como método de Adams-Bashforth de 2a ordem, e representa um método explícito onde é necessário conhecer a solução nos dois pontos anteriores. Usando polinômios de maior ordem, pode-se obter métodos de maior ordem. Um dos mais utilizados é o de quarta ordem: yn+1 = yn + ∆t 24 (55fn − 59fn−1 + 37fn−2 − 9fn−3) Neste caso, é necessário conhecer 4 pontos anteriores. Em comparação com outros métodos de quarta ordem, como RK4, o método de Adams-Bashforth costuma apresentar melhores resultados. 2.4.2. Método de Adams-Moulton O método de Adams-Moulton é uma modi�cação do método de Adams-Bashftorh que con- siste em avaliar o polinômio interpolador p(t) em tn+1, tn, tn−1, . . . ao invés de tn, tn−1, tn−2, . . .. 35 Por exemplo, considere novamente um polinômio de primeira ordem p(t) = αt+ β. Neste caso, as constantes são avaliadas através da relações: αtn+1 + β = fn+1 αtn + β = fn De onde se obtém que: α = fn+1 − fn ∆t β = fntn+1 − fn+1tn ∆t De modo que yn+1 pode ser obtido como: yn+1 = yn + ∆t 2 fn + ∆t 2 fn+1 Considerando que fn+1 = f(tn+1, yn+1), esta relação representa uma fórmula implícita de segunda ordem. Utilizando polinômios de ordem maior, obtém-se formulações de ordem superior. O mé- todo de Adams-Moulton de quarta ordem é dado como: yn+1 = yn + ∆t 24 (9fn+1 + 19fn − 5fn−1 + fn−2) Para evitar as di�culdades associadas ao uso de métodos implícitos, é comum adotar uma abordagem preditiva-corretiva, utilizando Adams-Bashforth como uma etapa preditiva e Adams-Moulton como corretora. Por exemplo, utilizando uma abordagem de segunda ordem, a etapa preditiva é dada como: y∗n+1 = yn + 3 2 ∆t · fn − 1 2 ∆t · fn−1 Sendo este valor posteriormente corrigido com a relação de Adams-Moulton: yn+1 = yn + ∆t 2 fn + ∆t 2 f ∗n+1 onde f ∗n+1 = f(tn+1, y ∗ n+1) Exemplo 03:) Utilizando uma abordagem preditiva-corretiva, utilize os métodos de Adams de segunda ordem para obter a solução aproximada do seguinte PVI até t = 0.4, utilizando ∆t = 0.1: dy dt = −y + 2t y(0) = 2 36 Esta equação foi avaliada anteriormente com o método de Runge-Kutta. Para utilizar Adams-Bashforth, neste caso é necessário conhecer a solução em dois pontos anteriores. Pode-se utilizar a condição inicial como um destes pontos: t0 = 0 y0 = 2 e o valor obtido com o método de Runge-Kutta para y1: t1 = 0.1 y1 = 1.82 Com isso, pode-se determinar o valor preditivo em t2 = 0.2: y∗2 = y1 + 3 2 ∆t · f(t1, y1)− 1 2 ∆t · f(t0, y0) = 1.577 Utilizando Adams-Moulton para corrigir o valor: y2 = y1 + ∆t 2 f(t1, y1) + ∆t 2 f(t2, y ∗ 2) = 1.68015 Repetindo o procedimento para o próximo passo de tempo (t3 = 0.3): y∗3 = y2 + 3 2 ∆t · f(t2, y2)− 1 2 ∆t · f(t1, y1) = 1.48313 y3 = y2 + ∆t 2 f(t2, y2) + ∆t 2 f(t3, y ∗ 3) = 1.57199 Finalmente, para o último passo de tempo (t4 = 0.4): y∗4 = y3 + 3 2 ∆t · f(t3, y3)− 1 2 ∆t · f(t2, y2) = 1.4902 y4 = y3 + ∆t 2 f(t3, y3) + ∆t 2 f(t4, y ∗ 4) = 1.4889 37 2.5. Estabilidade dos Métodos de Resolução de EDO's Na passagem do método de Euler explícito para métodos mais complexos, o principal objetivo foi obter métodos com melhor precisão (reduzir os erros de truncamento). No entanto, em muitos casos os resultados obtidos não só possuem uma baixa precisão como também são catastro�camente distintos da solução exata. Por exemplo, considere o seguinte PVI: dy dt = −12y y(0) = 1 A utilização do método de Euler explícito com ∆t = 0.2 gera o resultado ilustrado a seguir. Neste caso, o erro aumenta rapidamente conforme se avança na solução. Este problema está associado com a falta de estabilidade do método empregado. Antes de avaliar a estabilidade dos métodos, é necessário de�nir dois conceitos relacionados a análise de estabilidade: Consistência: Um método de passo único (como os da família de Runge-Kutta) é dito consistência se o erro de truncamento local ei tende a zero conforme ∆t → 0 para todos os passos de tempo, ou seja: lim ∆t→0 max 1≤i≤N |ei| = 0 onde N representa o número total de iterações. Convergência: Um método de passo único é dito convergente com respeito ao pro-blema de valor inicial dy dt = f(t, y) y(t0) = y0 se a diferença entre a solução obtida através do método numérico yi e a solução exata no mesmo ponto φi tender a zero conforme ∆t→ 0: lim ∆t→0 max 1≤i≤N |yi − φi| = 0 38 Dessa forma, um método consistente possui a propriedade de que as equações obtidas para cada passo se aproximada da própria equação diferencial conforme o passo de tempo tende a zero, de modo que o erro de truncamento local também tende a zero. Porém, além do erro de truncamento, deve-se considerar a in�uência do erro de arredonda- mento, associado ao fato de que os valores numéricos não são representados de forma exata. Na prática, os parâmetros do sistema, as condições inicias e todo o processo aritmético sub- sequente possuem erros associados com a precisão numérica �nita dos valores empregados. Para garantir que um determinado método seja convergente, deve-se então garantir que além de ser consistente, o erro de arredondamento deve �car limitado a um valor aceitável. O controle do erro de arredondamento está associado com a estabilidade do método. Este conceito é muito similar ao conceito de condicionamento de um sistema linear, no sentido de que um método é dito estável quando pequenas variações nos parâmetros ou condições iniciais levam a igualmente pequenas variações na solução obtida. A seguir será apresentada uma análise da estabilidade para os métodos de passo único. Para os métodos de passos- múltiplos, como existem diversas etapas de aproximação envolvidas em cada passo, a análise é mais complexa. 2.6. Análise de Estabilidade de Métodos de Passo Único Para avaliar a estabilidade de um dado método, considere o seguinte PVI: dy dt = f(t, y) y(t0) = y0 Considere que a solução exata deste PVI em um ponto tn é dada por φ(tn) e que a solução obtida com o uso de método numérico neste ponto é yn. Como comentado anteriormente, na resolução computacional do problema, existem associados erros de arredondamento devido à precisão limitada da representação numérica. Considere que y∗n é o valor arredondado obtido em uma máquina real e yn é o valor com precisão in�nitaobtido em uma máquina ideal. O erro total associado será a diferença entre a solução exata e o valor fornecido pelo computador: Erro total = φ(tn)− y∗n = (φ(tn)− yn) + (yn − y∗n) O primeiro termo representa o erro de truncamento e o segundo o erro de arredondamento associado ao ponto tn. Para que um dado método numérico seja adequado, são necessárias duas condições: 39 - O erro de truncamento acumulado (global) deve tender a zero conforme ∆t→ 0 (consis- tência); - O erro de arredondamento acumulado (que não pode ser eliminado) deve ser pequeno em comparação com a solução exata (estabilidade). 2.6.1. Equação Modelo: Decaimento de Primeira Ordem Para ilustrar as características de estabilidade de um método, considere a equação: dy dt = −λy y(0) = 1 onde λ > 0. Esta equação representa uma taxa de decaimento de primeira ordem, que surge, por exemplo, na análise da variação na fração de um dado reagente em uma reação de primeira ordem. Considerando novamente que φ(t) é a solução exata da equação, o valor calculado nume- ricamente pode ser expresso como: y(t) = φ(t) + ε(t) Substituindo na equação diferencial: d(φ+ ε) dt = −λ(φ+ ε) Como as solução exata deve satisfazer a equação diferencial, temos que: dε dt = −λε Esta equação serve como uma relação para determinar como o erro varia ao longo do tempo, ou seja, ao longo dos passos avaliados. Se o erro diminuir com o tempo, o método é estável, caso contrário será instável. Por exemplo, considere que o método de Euler explícito seja utilizado para avaliar ε(t): εn+1 = εn + ∆t(−λεn) = εn(1−∆tλ) → εn+1 εn = 1−∆tλ Para garantir a estabilidade, é preciso que o erro em tn+1 seja menor que o erro em tn, assim: ∣∣∣εn+1 εn ∣∣∣ ≤ 1 40 Considerando a relação anterior: |1−∆tλ| ≤ 1 Desse modo, o passo de tempo para garantir a estabilidade deve ser maior que zero e menor que 2/λ. A região contendo os valores de λ∆t que levam a uma solução estável é chamado de domínio de estabilidade linear. O parâmetro λ pode ser um número complexo, de modo que o domínio de estabilidade costuma ser representado no plano complexo. De�nindo z = −λ∆t, temos que para o método de Euler explícito ser estável a seguinte condição deve ser satisfeita: |1 + z| ≤ 1 Fazendo z = a+ bi: |1 + (a+ bi)| ≤ 1 → |(a+ 1) + bi| ≤ 1 Para um número complexo qualquer α + βi, o módulo é de�nido como: |α + βi| = √ α2 + β2 Assim, para o caso anterior, temos que: |(a+ 1) + bi| = √ (a+ 1)2 + b2 ≤ 1 → (a+ 1)2 + b2 ≤ 1 Esta relação representa um círculo de raio menor ou igual a 1, deslocado em uma unidade para a esquerda no eixo equivalente a a. Como a é a parte real de z e b a parte imaginária de z, o domínio de estabilidade do método de Euler explícito representa a região indicada na �gura a seguir. Considere agora que seja empregado o método de Euler implícito: εn+1 = εn + ∆t(−λεn+1) → εn+1(1 + ∆tλ) = εn Com isso: εn+1 εn = 1 1 + ∆tλ Como λ > 0 e ∆t > 0, esta relação mostra que o método de Euler implícito é estável para qualquer valor de ∆t, o que é uma característica comum dos métodos implícitos. 41 Considere agora a modi�cação de segunda ordem baseada na regra do ponto médio: εn+1 = εn + ∆t 2 (−λεn − λεn+1) De modo que: εn+1(1 + ∆tλ/2) = εn(1−∆tλ/2) → εn+1 εn = 1−∆tλ/2 1 + ∆tλ/2 Considerando a condição para estabilidade:∣∣∣1−∆tλ/2 1 + ∆tλ/2 ∣∣∣ ≤ 1 Neste caso, a solução também será estável para qualquer ∆t > 0. Para os demais métodos explícitos, a estabilidade também está condicionada a um domínio especí�co. Em particular para o caso dos métodos de RK de ordem maior que 1, pode-se mostrar que o domínio de estabilidade linear engloba o domínio associado ao método de Euler explícito (RK de primeira ordem). Assim, se a condição ∆t < 2/λ for respeitada, os métodos de RK de ordem superior serão estáveis. O domínio de estabilidade de métodos de Runge-Kutta de ordem 1 até 5 é representado na �gura a seguir: 2.7. Problemas Rígidos (Sti�) Algumas equações ou sistemas de equações diferenciais são classi�cados como rígidos (sti�). Não existe uma de�nição precisa para classi�car uma EDO como rígida, mas es- tes tipos de problema compartilham características em comum: 42 - Normalmente existem termos que levam à uma rápida variação na solução; - Problemas rígidos possuem uma variedade de escalas de tempo associadas, ou seja, em determinados pontos a solução varia muito mais rapidamente que em outras; - Métodos explícitos só são estáveis para a resolução de EDO's rígidas se o passo de tempo utilizado for muito pequeno; - Usualmente, o passo de tempo utilizado para garantir a estabilidade é menor que o necessário para garantir a convergência desejada. A equação teste utilizada anteriormente: dy dt = −λy y(t0) = y0 é um exemplo de equação rígida, especialmente para altos valores de λ. Esta equação possui solução da forma: y = y0e −λt 43 Como visto anteriormente, se um método explícito for empregado, deve-se usar um valor de ∆t su�cientemente pequeno para garantir a estabilidade da solução. Por exemplo, para o método de Euler explícito: ∆t < 2 λ Caso for necessário avaliar a solução até um tempo �nal t1, o número de passos necessários (n) será: n = t1 ∆t → n > λt1 2 Assim, o número de passos mínimo necessários é diretamente proporcional ao valor de λ. 2.7.1. Sistemas de Equações Diferenciais Rígidas Lineares Considere o seguinte PVI: dy1 dt = −500.5y1 + 499.5y2 y1(0) = 2 dy2 dt = 499.5y1 − 500.5y2 y2(0) = 1 A solução deste problema é a seguinte: y1(t) = 1.5e −t + 0.5e−1000t y2(t) = 1.5e −t − 0.5e−1000t Assim, a solução possui dois termos: um termo e−t que varia lentamente com o tempo e outro e−1000t que varia rapidamente. Para garantir a estabilidade, é preciso que a solução se mantenha estável durante toda a resolução, por isso, a estabilidade deve ser assegurada para o termo e−1000t. De modo geral, os valores de ∆t empregados devem ser determinados sem o conhecimento da solução. O PVI anterior pode ser reescrito como: dY dt = −500.5 499.5 499.5 −500.5 Y = AY O valor de ∆t necessário é de�nido com base nos autovalores da matriz A. Lembrando que os autovalores λ são de�nidos como os valores onde: |A− λI| = 0 44 Assim, para este caso:∣∣∣∣∣−500.5− λ 499.5499.5 −500.5− λ ∣∣∣∣∣ = 0 → (−500.5− λ)2 − 499.52 = 0 Resolvendo a equação de segundo grau, obtém-se as seguintes raízes: λ1 = −1 λ2 = −1000 Para garantir a estabilidade em um método explícito, pode-se considerar a função teste usada anteriormente dy/dt = −λy, com λ sendo o maior dos autovalores (em módulo) associados com o problema. Por exemplo, caso o método de Euler explícito seja empregado, deve-se garantir que: ∆t ≤ 2 |λmax| Considerando que no exemplo |λmax| = 1000, então temos que o passo de tempo máximo é de 0.002. Caso o menor autovalor fosse empregado, seria obtido um valor de ∆tmax = 2. Após os instantes iniciais, o problema passa a ser controlado pelo menor autovalor. Para medir a importância da rigidez na resolução do problema, pode-se determinar o grau de rigidez (sti�ness ratio) do problema, de�nido com: SR = |λmax| |λmin| Normalmente, para SR > 20 o problema já é classi�cado como rígido, sendo necessário avaliar com cuidado os passos de tempo empregados. 2.7.2. Sistemas Não-Lineares Para problemas não-lineares, a análise é mais complexa. Para problemas autônomos, da forma: dY dt = f(Y ) pode-se linearizar a equação através de uma expansão em série de Taylor em torno do ponto tn. Desconsiderando os termos de alta ordem: dY dt = f(Yn) + J(tn)(Y − Yn) onde J(tn) é a matriz Jacobiana avaliada em t = tn, de�nida como: J(tn) = aij = [ ∂fi(y) ∂yj ] tn 45 Por exemplo, em um sistema 2× 2 da forma: dy dt = f(x, y) dx dt = g(x, y) O Jacobiano em um ponto tn é determinado da forma: J(tn) = ∂f∂x ∣∣∣ tn ∂f ∂y ∣∣∣ tn ∂g ∂x ∣∣∣ tn ∂g ∂y ∣∣∣ tn Neste caso, a razão de rigidez é de�nida com base nos autovalores da matriz Jacobiana. Como J(tn) pode depender do tempo, a razão de rigidez podevariar ao longo da solução. Os valores de ∆t necessários para garantir a estabilidade também podem ser de�nidos com base nos autovalores da matriz Jacobiana. Obs.: Caso tn for um ponto �xo, os autovalores do Jacobiano também servem para de�nir se o ponto é estável ou não, sendo instável caso algum autovalor tenha parte real positiva. Exemplo: Considere o seguinte PVI: dy dt = y2x y(0) = 2 dx dt = −400yx2 x(0) = 1 Determine a razão de rigidez em t = 0. Caso o método de Euler explícito for empregado, qual o valor máximo de ∆t? 46 3. Métodos Numéricos para Problemas de Valor de Contorno Equações diferenciais de ordem maior que um podem gerar problemas de valor inicial (PVI) ou problemas de valor de contorno (PVC), dependendo da forma como as condições conhecidas são especi�cadas. Por exemplo, considere a EDO: y′′ + p(x)y′ + q(x)y = g(x) Até o momento, foram analisados principalmente casos onde condições iniciais conhecidas são especi�cadas da forma: y(x0) = y0 y ′(x0) = y ′ 0 originando desta forma um PVI. Para a resolução numérica de PVI's, pode-se partir da condição inicial e ir avançando até um tempo �nal arbitrário. Em muitos casos, os problemas envolvem condições conhecidas em pontos diferentes, sendo estas chamadas de condições de contorno, podendo ser expressas, por exemplo, como: y(α) = y0 y(β) = y1 De forma geral, os PVC's envolvem uma coordenada espacial como variável independente. Assim, a resolução de um PVC's consiste em buscar uma solução que satisfaz a equação diferencial no intervalo α < x < β juntamente com as condições de contorno especi�cadas. Isto implica que existem duas condições, em pontos diferentes do domínio de solução, que deve ser simultaneamente satisfeitas. Por isso, os métodos de marcha (como os de Runge- Kutta) não podem ser empregados neste caso. 47 3.1. Estratégias de Solução de PVC's Os métodos para resolver problemas de valor de contorno se dividem em duas categorias: os baseados em transformar o PVC em um PVI e os baseados em discretizar a equação utilizando métodos de diferenças �nitas. Os métodos que transformação PVC's e PVI's consistem basicamente em utilizar uma das condições de contorno como condição inicial e assumir (chutar) diferentes valores para uma segunda condição inicial de modo que o resultado obtido satisfaça a segunda condição de contorno. Por exemplo, considere a seguinte equação utilizada para descrever a variação na temperatura em uma barra que perde calor para o ambiente por convecção, como apresentado na �gura a seguir. Considerando que a barra seja muito �na, com raio muito menor que o comprimento, pode-se assumir que a equação que descreve a variação na temperatura ao longo de x pode ser expressa como: d2T dx2 + h′(Ta − T ) = 0 onde h′ é um coe�ciente de troca térmica e Ta é a temperatura ambiente. As condições de contorno corresponde a temperatura �xas nas extremidades, em x = 0 e x = L, de modo que: T (0) = T1 T (L) = T2 Para transforma a equação em um PVI equivalente, primeiramente é preciso aplicar o método de redução de ordem. Escrevendo o PVI como: dT dx = z T (0) = T1 dz dx + h′(Ta − T ) = 0 z(0) = z0 48 O valor de z0 não é conhecido, pois o valor da derivada de T em x = 0 não foi especi�cado. O método utilizado neste caso consiste em assumir diferentes valores para z0 até que a con- dição T (L) = T2 seja satisfeita. Para equação lineares, este método funciona razoavelmente bem, pois pode-se interpolar os valores para T (L) obtidos para diferentes valores de z0 para encontrar o valor de z0 que corresponde a solução correta do problema. No entanto, de modo geral este método é pouco utilizado, em particular porque os métodos baseados em aproximações por diferenças �nitas costumam ser mais e�cientes, como ser apresentado a seguir. 3.2. Aproximações por Diferenças Finitas O método de diferenças �nitas é um método de discretização de equações diferenci- ais. Isto signi�ca que ele transforma uma função contínua em uma representação dis- creta (pontos). Por exemplo, considere uma função f(x) = 2x de�nida em um intervalo entre 0 e 1. Esta função pode ser representada de forma discreta como, por exemplo, f [x] = [0, 0.4, 0.8, 1.2, 1.6, 2], assumindo um espaçamento entre os pontos de 0.2, ou seja, esta representação está relacionada como um domínio discreto da forma x = [0, 0.2, 0.4, 0.6, 0.8, 1]. As soluções obtidas com a aplicação do método de diferenças �nitas sempre serão discretas. Caso for necessário obter os valores para algum valor de x que não corresponda exatamente aos pontos do domínio, pode-se interpolar os valores. A primeira etapa da aplicação do método de diferenças �nitas consiste exatamente em de�nir o domínio discreto onde a solução será buscada. Por exemplo, considere o caso apresentado anteriormente para a distribuição de calor em uma barra estacionária. A região onde se deseja obter a distribuição de temperatura é no intervalo de x = 0 até x = L. Este intervalo corresponde ao domínio de solução da equação diferencial. No entanto, ele está em uma forma contínua e não discreta. Para discretizar o domínio, deve-se dividi-lo em um determinado número de pontos. Por exemplo, considere que L = 1 e que se deseja dividir o domínio em pontos com espaçamento ∆x = 0.1, ou seja, deseja-se dividir o domínio em 10 elementos de igual tamanho. Quanto mais elementos forem utilizados, maior será a precisão do método, porém o gasto computacional também irá aumentar. O domínio discreto será então: x = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9.1] 49 Como pode ser visto, este conjunto contém 11 pontos. De maneira geral, o número de pontos sempre será igual ao número de elementos em que o domínio é dividido mais 1. Este vetor pode ser representado de uma forma mais simples como: x[i] = i∆x i = 0, 1, 2, . . . 10 Na resolução de PVI's, não é necessário inicialmente de�nir o domínio de solução, pois pode-se continuar avançando por quantos passos forem necessários, partindo de um valor inicial. No entanto, para o caso de PVC's, o domínio de solução é fechado e deve ser considerado como um todo. Dependendo da formulação utilizada, pode até mesmo ser necessário obter os valores para todos os pontos ao mesmo tempo. A estratégia do método de diferenças �ntas consiste em buscar equações algébricas que aproximem a solução em cada ponto i. Para isso, as derivadas são aproximadas como relações algébricas envolvendo a solução em diferentes valores de i. Esta aproximação pode ser realizada de diferentes formas, dependendo da precisão desejada e da natureza do problema e das condições de contorno. Porém, a origem destas aproximações sempre é uma aproximação em série de Taylor em torno de cada ponto i, como será discutido a seguir. 3.2.1. Aproximação da Derivada Primeira Para apresentar o processo de discretização da derivada de uma função contínua y(x) em um intervalo 0 ≤ x ≤ L, será considerado um domínio discreto como o apresentado na �gura a seguir, onde o domínio físico contínuo (região entre 0 e L) é dividido em N + 1 pontos. Lembrando novamente, o objetivo do método de diferenças �nitas é obter aproximações para o valor da função y(x) em cada um destes N + 1 pontos. Esta representação discreta do domínio de solução é normalmente chamada de grid nu- mérico ou malha numérica. Como visto em aulas anteriores, a expansão em série de Taylor pode ser utilizada para avaliar o valor de uma função em um dado ponto com base no valor conhecido em outro ponto. Dependendo da forma como a aproximação é realizada, obtém-se diferentes formulações para o método de diferenças �nitas, como será apresentado a seguir. 50 Aproximação para Frente (Foward) Considere que se deseja aproximar o valor em xi+1 com base no valor em xi, ou seja, deseja-se aproximar y(xi+1) = yi+1 com base em y(xi) = yi. A expansão em série de Taylor neste caso pode ser expressa como: yi+1 = yi + (xi+1 − xi) dy dx ∣∣∣ x=xi + (xi+1 − xi)2 2! d2y dx2 ∣∣∣ x=xi+ (xi+1 − xi)3 3! d3y dx3 ∣∣∣ x=xi + . . . (3.1) Considerando novamente que (xi+1 − xi) seja relativamente pequeno, os termos de alta proporcionais a (xi+1 − xi)2, (xi+1 − xi)3, . . . podem ser desprezados. Assim, a derivada primeira da função y(x) no ponto xi pode ser aproximada como: dy dx ∣∣∣ x=xi = yi+1 − yi xi+1 − xi De�nindo ∆x = xi+1 − xi, a expressão pode ser dada por: dy dx ∣∣∣ x=xi = yi+1 − yi ∆x Esta expressão é conhecida como aproximação por diferenças �nitas para frente (ou método de diferenças �nitas para frente), pois utiliza o valor da função em um ponto a frente xi+1 para estimar o valor da derivada em um ponto anterior xi. Fazendo o limite de ∆x tendendo a zero: dy dx ∣∣∣ x=xi = lim ∆x→0 yi+1 − yi ∆x obtém-se a própria de�nição da derivada em um ponto. Portanto, esta formulação é con- sistente. Como pode ser visto, este método é equivalente ao método de Euler explícito, sendo que a diferença é que no método de Euler deseja-se aproximar o valor da função em diferentes pontos com base no conhecimento da derivada dy/dt = f(t, y) e no método de diferenças �nitas deseja-se aproximar o valor da derivada com base no valor em diferentes pontos. Nos dois casos, porém, ocorre uma linearização da função em torno de um ponto. Na aproximação da derivada, foram desprezados os termos O(xi+1− xi)2, assim, o erro de truncamento local do método de diferenças para frente é da ordem de O(∆x2). De forma semelhante ao apresentado para o método de Euler, pode-se mostrar que quando aplicado para avaliar N pontos, o erro associado será da ordem de O(xi+1 − xi) = O(∆x), ou seja, a aproximação por diferenças para frente é um método de primeira ordem. 51 Aproximação para Trás (Backward) De forma semelhante ao realizado para obter a derivada em xi com base na expansão para obter yi+1, pode-se realizar uma expansão para obter a função em xi−1: yi−1 = yi + (xi−1 − xi) dy dx ∣∣∣ x=xi + (xi−1 − xi)2 2! d2y dx2 ∣∣∣ x=xi + (xi−1 − xi)3 3! d3y dx3 ∣∣∣ x=xi + . . . Considerando que o espaçamento entre os pontos é constante, temos novamente que ∆x = xi − xi−1 = −(xi−1 − xi), assim: yi−1 = yi − (∆x) dy dx ∣∣∣ x=xi + (∆x)2 2! d2y dx2 ∣∣∣ x=xi − (∆x) 3 3! d3y dx3 ∣∣∣ x=xi + . . . (3.2) Novamente, desprezando os termos de ordem maior ou igual a 2, obtém-se: dy dx ∣∣∣ x=xi = yi − yi−1 ∆x esta aproximação é conhecida como aproximação por diferenças �nitas para trás (ou mé- todo de diferenças �nitas para trás), e também representa uma aproximação de primeira ordem para a derivada em um dado ponto xi. A utilização dos métodos para trás ou para frente é, a princípio, equivalente. A única restrição ocorre nos pontos extremos do domínio. Por exemplo, no ponto x0 não pode ser aplicado o método para trás pois não existe um ponto anterior a este. De forma semelhante, no ponto xN a formulação para frente não pode ser utilizada, pois de maneira equivalente não existe nenhum ponto após este para ser utilizado de base. Aproximação Central De forma geral, a utilização de métodos de primeira ordem para a discretização de todos os pontos do domínio não costuma apresentar bons resultados, especialmente nos casos envolvendo gradientes aproximadamente simétricos em relação a direção x, como por exemplo problemas envolvendo condução de calor ou difusão de massa. Uma aproximação de segunda ordem pode ser obtida fazendo a expansão para yi+1 (Eq. 1) menos a expansão para yi−1 (Eq. 2). Com isso, obtém-se: yi+1 − yi−1 = 2∆x dy dx ∣∣∣ x=xi + 2(∆x)3 3! d3y dx3 ∣∣∣ x=xi + . . . Desprezando agora os termos da ordem de O(∆x)3, pode-se obter a seguinte expressão para uma aproximação da derivada primeira em xi: dy dx ∣∣∣ x=xi = yi+1 − yi−1 2∆x 52 Esta expressão é conhecida como aproximação pode diferenças �nitas central (ou método de diferenças �nitas central). Neste caso, o erro de truncamento local é da ordem de O(∆x)3, de modo que o erro global será da ordem de O(∆x)2. Assim, esta aproximação é de segunda ordem. Existem aproximações de ordem superior, porém a aproximação central é a mais empre- gada, sendo adequada para a maioria dos casos. Novamente, dever-se observar que esta aproximação não pode ser aplicada nos pontos extremos do sistema. A melhor estratégia para a discretização da equação é utilizar o método central para os pontos internos, o método para frente no ponto x0 e o método para trás no ponto xN . Uma comparação geométrica dos três esquemas de discretização é apresentada na �gura a seguir. 3.2.2. Aproximação da Derivada Segunda O método de diferenças �nitas pode ser aplicado também para a discretização de derivadas de maior ordem. As aproximações para a derivada segunda são obtidas considerando a de�nição da derivada em um ponto. Da mesma forma que a derivada primeira pode ser aproximada em termos da variação da função em dois pontos, a derivada segunda pode ser aproximada em termos da variação na derivada primeira em dois pontos. Por exemplo, utilizando um esquema para frente: d2y dx2 ∣∣∣ x=xi = dy dx ∣∣∣ x=xi+1 − dy dx ∣∣∣ x=xi xi+1 − xi 53 onde a derivada avaliada no ponto xi+1, utilizando novamente um esquema para frente, é obtida de forma equivalente a derivada no ponto xi como: dy dx ∣∣∣ x=xi+1 = yi+2 − yi+1 xi+2 − xi+1 Considerando que os pontos sejam igualmente espaçados (grid igualmente espaçado), pode-se juntar as expressões para de derivada em xi+1 e a expressão para a derivada em xi, obtém-se a seguinte expressão para a derivada segunda em xi: d2y dx2 ∣∣∣ x=xi = yi − 2yi+1 + yi+2 (∆x)2 Esta expressão representa o esquema para frente aplicado à derivada segunda. Como ele se baseia em esquemas de primeira ordem, também é uma aproximação de primeira ordem. Fazendo um procedimento semelhante utilizando o esquema para trás, obtém-se: d2y dx2 ∣∣∣ x=xi = yi − 2yi−1 + yi−2 (∆x)2 sendo este o esquema para trás aplicado para a derivada segunda, sendo também um método de primeira ordem. Utilizando o esquema central, obtém-se a seguinte aproximação: d2y dx2 ∣∣∣ x=xi = yi+1 − 2yi + yi−1 (∆x)2 sendo esta uma aproximação de segunda ordem. Novamente, a estratégia que normalmente resulta em um melhor resultado é utilizar o esquema central para os pontos internos e os esquemas para frente e para trás para o primeiro e para o último ponto, respectivamente. 3.2.3. Discretização das Condições de Contorno Além da equação diferencial envolvendo a derivada da função, a resolução de Problemas de Valor de Contorno necessita que determinadas condições de contorno também sejam em determinados pontos do domínio. Para o caso de condições que consistem em de�nir o valor da função em um dado ponto, basta atribuir este valor para a função discretizada no ponto. Por exemplo, considere o exemplo anterior de transferência de calor em uma barra, onde a temperatura nas duas extremidades era conhecida: T = T1 em x = 0 e T = T2 em x = L. Considerando que 54 x0 corresponde ao ponto inicial x = 0 e xN corresponde ao ponto x = L, as condições de contorno são implementadas de�nindo: T0 = T1 TN = T2 Em muitos casos, a condição de contorno envolve a própria derivada de função em algum ponto. Por exemplo, considere que no exemplo da transferência de calor na barra metálica, a extremidade x = L fosse mantida isolada termicamente. Neste caso, a condição de contorno é expressa como: dT dx ∣∣∣ x=L = 0 Neste caso, é necessário discretizar a condição de contorno para transforma-la numa relação algébrica. Como este ponto corresponde ao último ponto do domínio, os esquemas para frente e central não pode ser utilizados, pois iriam depender de um ponto yN+1 que está fora do domínio de solução. Assim, é necessário utilizar o esquema para trás, de modo que a condição pode ser discretizada como: dT dx ∣∣∣ x=L = TN − TN−1 xN − xN−1 = 0 → TN = TN−1 ou seja, a condição de derivada nula na posição x = L implica que a variável em no ponto xN é igual a variável no ponto anterior
Compartilhar