24_-_CalcNum__-_EDO_s_-_métodos_numéricos.pdf

Transcript

1. LNCC Métodos numéricos para EDO’s Prof. Paulo R. G. Bordoni UFRJ
2. Nas transparências a seguir, mostrarei que existem diversas possibilidades de obter aproximações numéricas para soluções de uma EDO.
3. Na realidade obteremos aproximações para a solução do problema de valor inicial ቊ ′ = (, ) 0 = 0 . Vamos assumir que a solução ↦ () do PVI é única e está definida num intervalo num intervalo [, ].
4. Vamos criar uma partição , = = 0 , 1 , … , = e calcular aproximações 0 , 1 , … , para os valores exatos 0 , 1 , … , ( ) da solução. 1 0 0 1 1
5. Começamos observando que a condição inicial 0 = 0 fornece o ponto de partida da solução ↦ (). Sim Mestra, mas só isso. A solução () é como o rabo de um cachorro, preso no ponto 0 , 0 − abana prá cima e prá baixo! 0 0 () (0 , 0 )
6. 0 0 () ′(0 ) É mesmo! Isto fixa a direção do rabo do cachorro (da reta tangente à curva), como nos perdigueiros quando apontam uma codorna na moita! Não apenas, Surfista! A equação ′ = (, ) fornece a inclinação inicial da curva ↦ (), visto que ′ 0 = (0 , 0 ) . Em outras palavras, seu cachorro para de abanar o rabo
7. 0 ′(0 ) 1 1 = (1 ) 0 Portanto a equação da reta tangente ↦ (), com inclinação ′ 0 e passando pelo ponto (0 , 0 ), é = 0 + ′(0 )( − 0 ). Em particular calculamos 1 = 1 por 1 = 0 + ′(0 )(1 − 0 ). reta tangente
8. 0 ′(0 ) 1 1 = (1 ) 0 E como a reta tangente ↦ () “cola” na curva pertinho do ponto de tangência poderemos aproximar 1 − na curva ↦ − pelo valor 1 = (1 ) Sim Cabelos de Fogo, mas estaremos usando a ideia de colar ao contrário. Agora é a curva que “cola” na reta tangente! Aproximação (1 )
9. 0 ′(0 ) 1 1 0 ℎ Portanto sabemos como calcular 1 = 0 + ′(0 )(1 − 0 ) e assim a aproximar a solução ↦ () do problema de valor inicial num ponto 1 = 0 + ℎ quando ℎ é bem pequeno. Basta colocar 1 ≅ 1 = 0 + ′(0 ) ∙ ℎ 1 ≅ 1 (1 )
10. Mestra, pela figura parece que o erro é muito grande! Nem tanto Loirinha, comparando 1 com a série de Taylor da solução () para = 0 + ℎ, obtemos 0 + ℎ = 0 + ′ 0 ℎ + ′′ 0 ℎ2 2 + … d’onde concluímos que o erro nessa aproximação é ℴ(ℎ2).
11. Sim, eu e o Cabelos de a Fogo, já falamos que reta tangente e a curva “colam” uma na outra perto do ponto de tangência 0 , 0 . Isso mesmo Surfista. Para ℎ = 0.1 o erro é da ordem de 0.01, dez vezes menor.
12. E como continuamos, Mestra ? Bem, é só repetir a ideia trocando (0 , 0 ) por (1 , 1 ) e assim por diante. 0 ′(1 ) ≅ (1 , 1 ) 2 2 (2 ) 0 ℎ 1 1 (1 )
13. É. Mas o erro vai acumulando em cada passo. Sim, temos o erro local ℴ(ℎ2) e o erro acumulado − global, ℴ(ℎ).
14. Essas ideias devem ter sido a inspiração de L. Euler (circa 1768) para aproximar a solução do PVI nos pontos da partição , = = 0 , 1 , … , = : ቊ +1 = + ℎ ∙ , , = 0,1, ⋯ , − 1 0 = (0 )
15. No programa a seguir codifico o método de Euler para aproximar o PVI ቊ ′ = (, ) 0 = 0 numa partição , = = 0 , 1 , … , = de um intervalo [, ].
16. Marquei a entrada de dados do programa do Mestre:
17. E aqui mostrei a implementação do método de Euler:
18. Esse é o gráfico da solução do programa para o método de Euler comparado com o resultado da odeint( ) da Scipy:
19. Observamos que a fórmula de Euler utiliza as diferenças finitas avançadas. E isto nos leva a considerar as diferenças finitas retardadas: ( ), = ′( ) ≅ ( ) − (−1 ) ℎ
20. E nesse caso teremos a aproximação ≅ −1 + ℎ( , ) Ou, trocando por + 1: +1 ≅ + ℎ(+1 , +1
21. Então usando ainda as aproximações = e +1 = +1 teremos o chamado método retardado (backward) de Euler ቊ +1 = + ℎ ∙ +1 , +1 , = 0,1, ⋯ , − 1 0 = (0 ) Note Surfista, que esse esquema não fornece uma fórmula explícita para calcular +1 . Para utilizá-la precisamos resolver uma equação para obter +1 em termos de e isto depende inteiramente da natureza da função (, ).
22. Para o PVI do último exemplo onde, ′ = ∙ (), a equação a resolver é +1 = + ℎ ∙ +1 ∙ (+1 ) Ou +1 = 1 − ℎ( + ℎ) . Porém, para outras equações pode não ser tão simples explicitar +1 em termos de .
23. Vamos trabalhar uma outra forma de obter o método de Euler.
24. ′ ′() ′ Mestre, se é o tempo e () a posição de uma partícula no instante então ’() é a velocidade da partícula nesse instante. velocidade () Os gráfico da posição e velocidade em em função do tempo A associação entre a velocidade, posição e o tempo é descrita pela EDO: ′ = (, ) posição
25. Consideremos novamente uma partição uniforme , = = + ℎ para = 0, ⋯ , Então, integrando a equação ′ = (, ) de à +1 obtemos ׬ +1 ′ = ׬ +1 , ou +1 − = ׬ +1 , .
26. A integral do lado direito pode ser aproximada, grosseiramente, pela área azul indicada na figura, isto é ׬ +1 , () ≅ ( , ) × (+1 − ) Então: +1 − ≅ ( , ) ∙ ℎ ′ +1 ′( ) ′
27. ( , ) +1 (( ), ) ′ Pois é Mestra, para = 0 pela condição inicial 0 = 0 essa expressão fornece 1 ≅ 0 + (0 , 0 ) ∙ ℎ Exatamente a aproximação do método de Euler para (1 ) e que forneceu inspiração para o Euler deduzi-lo. Não vejo porquê repetir tudo novamente !
28. Bem, porque ela sugere outras aproxima- ções para a integral ׬ +1 , . Por exemplo pela área do trapézio azul: ( , ) +1 (( ), ) ((+1 ), +1 ) ′ Que fornecerá a fórmula aperfeiçoada de Euler, também conhecida como fórmula de Heun: +1 = + ℎ , + +1 , +1 2
29. Mas Mestre, com isto recairemos no mesmo problema da fórmula por diferenças retardadas. Possivelmente equações difíceis de resolver! Sim, mas podemos aproximar o valor de +1 em +1 , +1 pela fórmula de Euler: +1 ∗ = + ℎ ∙ ( , )
30. Com isso, aproximamos a fórmula de Heun em dois passos: i. 1º usamos a fórmula de Euler para aproximar +1 por +1 ∗ = + ℎ ∙ ( , ) i. Depois calculamos +1 = + ℎ , + +ℎ, +1 ∗ 2 usando a a inclinação média mostrada acima. ′ +1 ′ ′ ′ +1 ′ + ′ +1 2
31. (1 ) 0 0 1 1 ∗ ∡(1 , 1 ∗) ∡(0 0 ) 1 () Uma interpretação geométrica do método de Heun, através das inclinações (os ∡). ∡ médio = ′ +′ +1 2
32. A fórmula do RK é a seguinte: +1 = + ℎ + 2 + 2 + 6 onde = ( , ) = ( + ℎ 2 , + ℎ 2 ) = ( + ℎ 2 , + ℎ 2 ) = ( + ℎ , +ℎ ) O método de Runge-Kutta é uma média ponderada de inclinações nos extremos esquerdo, direito e o ponto médio , +1 e +1/2 do intervalo. Inspiração no método de Heun?
33. O RK é famoso e eficiente, seu erro local é (ℎ5) e o global ℎ4 , por isso chamado RK4. Loirinha, note que se = () então teremos = ( ), = + Τ ℎ 2 , = ( + Τ ℎ 2) e = ( + ℎ) e então +1 = ℎ 6 + 4 + Τ ℎ 2 + + ℎ , é a fórmula de Simpson para aproximar ׬ +1
34. Marquei tanto o pedido da solução exata e da EDO como o da entrada de dados A implementação do Runge- Kutta, da solução exata e da scipy.odeint( )
35. Aqui marquei o cálculo do RK 4:
36. Aqui desenhamos os três gráficos:
37. Os dados e a tabela comparativa:
38. Os gráficos da solução exata, do Runge-Kutta e da odeint( )
39. O que há na scipy sobre resolução de EDOs:
40. São apenas dois métodos, o odeint e o ode. Mas, conforme veremos, muito poderosos.
47. Surfista, considere um pêndulo de massa preso na ponta de uma haste de comprimento sob ação da gravidade , como na figura. Seu problema é descrever a variação temporal do ângulo = () do pêndulo (da haste) com a vertical. =
48. = O deslocamento (comprimento do arco) do pêndulo é = e a força tangencial que o faz retornar à posição de equilíbrio vertical é = − (). Logo − = = ′′ = ′′, pela 2ª lei de Newton. Portanto o movimento é descrito pela EDO não-linear de 2ª ordem ′′ = − Τ () Para tornar a coisa mais real vamos acrescentar um fator amortecimento à velocidade angular ′, de forma que a equação do movimento fica: ′′ + ′ + Τ () = 0.
49. O problema de valor inicial de 2ª ordem associado à oscilação do pêndulo fica ( → ) ൞ ′′ + ′ + Τ = 0, , > 0 0 = 0 ′ 0 = 0 ′ θ → porque não existe a letra grega em ASCII. Ele é equivalente ao sistema de 1ª ordem ′ ′ = − − Τ () (0) (0) = 0 0
50. Utilizamos a . ( ) para obter a solução. O programa está a seguir.
51. O programa. Marquei a entrada de dados.
52. Marquei a função que define o sistema de eqs. lineares e a chamada da odeint( )
53. Já o scipy.integrate.ode é mais que um método, trata-se de uma classe. Não vamos explorá-lo porque foi lançada uma nova API para equações diferenciais na scipy. Apenas vou mostrar o que há nele.
54. A nova API lançada em maio/2019
55. A classe scipy.integrate.ode
64. Tchau, esta foi a última aula