Pro Yearly is on sale from $80 to $50! »

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

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

4acc58a03aa964e2f04b538836f2d468?s=128

Paulo Bordoni

July 09, 2019
Tweet

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.
  41. None
  42. None
  43. None
  44. None
  45. None
  46. None
  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

  56. None
  57. None
  58. None
  59. None
  60. None
  61. None
  62. None
  63. None
  64. Tchau, esta foi a última aula