EDO's, Problemas de valor inicial

...

Transcript

1. UFRJ EDO’s, Problemas de valor inicial Prof. Paulo R. G. Bordoni
2. Vamos estudar agora métodos numéricos para resolver equações diferenciais ordinárias.
3. Vamos analisar cada componente dessa expressão em detalhe. Equações diferenciais ordinárias de são descritas, taquigraficamente, pela expressão ′ = (, ).
4. Se = 1 2 ⋯ e = 1 2 ⋯ são vetores de ℝ então ′ = (, ) é uma equação vetorial que descreve um sistema de EDOs. O estudo de EDO’s remonta aos inventores do Cálculo Diferencial e Integral, Newton e Leibnitz, entorno dos anos 1670. Euler (século XVIII) apresentou contribuições fundamentais.
5. Na forma de equações o sistema fica ′1 = 1 (, ) ′2 = 2 (, ) ⋮ ′ = (, ) , com = 1 2 ⋮ . Sistema de EDOs, Mestra? Francamente, não entendi.
6. Sim, basta escrevermos ቊ ′ = ′ = (, , ) Uma equação diferencial ordinária de 2ª ordem, ′′ = (, ′, ), se reduz a um sistema de duas equações de 1ª ordem.
7. Por exemplo:
8. Se a EDO não depende explicitamente de , isto é, ′ = () o sistema é dito autônomo. Já no caso da ser independente de , ou seja ′ = (), pelo teorema fundamental do cálculo, a EDO se reduz a uma integral indefinida: − (0 ) = ׬ 0 . Se , = + (), com , funções de , a EDO ′ = (, ) é linear. Caso = 0 a equação é homogênea. Para sistemas, é uma matriz × e um vetor.
9. Mais exemplos:
10. O primeiro dos ingredientes é uma função (, ) ↦ (, ), com um domínio ⊂ ℝ2 e a valores reais: (, ) ∈ ℝ. Físicos, engenheiros e meteorologistas usualmente chamam tais funções de campos escalares.
11. O CPTEC/IMPE usa curvas de nível para mostrar valores médios das temperaturas máximas e mínimas do campo de temperaturas no Brasil e regiões. Vejam:
12. Temperatura, a pressão, radiação solar, a pluviometria em cada ponto (do Brasil) são grandezas escalares, mensuradas pelo CPTEC/INPE. O CPTEC usa instrumentos sofisticados. Clique onde marquei para saber mais.
13. O gráfico de uma função : ℝ2 → ℝ e suas curvas de nível.: Esse exemplo é clássico!
14. Um outro exemplo. Agora uma sela.:
15. O código dos exemplos.
16. A continuação do código.
17. Apenas em si, (, ) define um campo escalar (, ) ↦ (, ) num domínio ⊂ ℝ2 e a valores reais. Entretanto, a igualdade em ′ = (, ) informa que o valor numérico desse campo (, ) ↦ (, ) deve ser interpretado como uma derivada ′ = Τ .
18. () ′ = (, ) 0 1 (1 , 1 ) 2 1 0 2 (0 , 0 ) (2 , 2 ) Se entendi, Galileu, o campo escalar (, ) ↦ (, ) define a inclinação ′(t) da função ↦ (), para cada ponto .
19. 0 1 (1 , 1 ) 2 1 0 2 (0 , 0 ) (2 , 2 ) Pois é Cabelos de Fogo, poderemos também pensar em um campo de vetores, se levarmos em consideração o sentido de percurso da curva ↦ (, ()). Repetindo a fala do Surfista, o campo (, ) ↦ (, ) é um campo escalar de inclinações definido em cada ponto (, ) ∈ ⊂ ℝ2.
20. Em outras palavras para cada ponto (, ) ∈ o campo associa um vetor , = (, ) (, ) . Imagine o fluxo d’água num rio e pense no campo de velocidades (das partículas de água). Um campo de vetores com domínio ⊂ ℝ2 é uma aplicação : → ℝ2 dada por (, ) ↦ , = (, ) (, )
21. Pois é, um campo escalar. Entretanto a igualdade ′ = (, ) força a entender os valores (, ) como inclinações, já que ′ = Τ é uma derivada. Mas a função (, ) ↦ (, ) assume valores reais, isto é, (, ) ∈ ℝ. Francamente, não entendo!
22. Este programa traça o campo de vetores definido pela função (, ) ↦ (, ). O rabinho deles fica preso nos pontos da meshgrid().
23. Os vetores são traçados com a função quiver() da matplotlib.
24. A função quiver() da matplotlib.
25. Uma execução do programa, com o campo definido por (, ) ↦ (, )
26. Modifiquei o programa, acrescentando o trecho recortado, onde os vetores do campo são normalizados com uma gradação de cores para definir seus tamanhos (normas). As cores mais escuras indicam vetores com norma maior .
27. Prefiro assim, Mestra, há mais informação visual.
28. Mais um exemplo:
29. 1 0 1 0 () () () Imagine Loirinha, que um ponto = () está se deslocando sobre um caminho , parametrizado por ↦ = [ , ], para ∈ [0 , 1 ] () é o vetor deslocamento, a posição em cada instante de tempo .
30. Para concretizar, Loirinha: você se deslocando de bicicleta na pista da Lagoa Rodrigo de Freitas. A Lagoa é a curva e sua posição no instante é = () () . O deslocamento do Surfista é outro!
31. 1 0 1 0 () () ′() () A velocidade de deslocamento do ponto () é dada por ↦ ′ = [ ′(), ′ ] É o vetor tangente à curva no ponto (). Eu ando mais rápido!
32. ′() ′() ′() () () () Claro Mestre, é a velocidade da minha bicicleta! Atenção Surfista: O vetor velocidade ′() e suas coordenadas ′() e ′ são referidas à um sistema de referência “amarrado” em ().
33. Fiz um programa que recebe a expressão da função (, ), calcula a solução ↦ () do PVI ቊ ′ = (, ) 0 = 0 e traça alguns vetores do campo (, ) ↦ (, ) tangentes ao gráfico de .
34. Os dados utilizados nessa execução do programa são:
35. Eis o código do programa. Aqui só mostro a entrada de dados
36. Usei a função odeint() da scipy.integrate para calcular a solução = ().
37. Usei uma escala de cores (crescendo do + claro ao + escuro) para pintar os vetores tangentes desenhados com a função arrow() da matplotlib.
38. A função odeint() da scipy.integrate.
39. 1 0 0 () 1 = [ , ] ′ = [ 1, ′ ] Para funções, a parametrização é ↦ [, ]. E portanto o vetor velocidade é ↦ ′ = [ 1, ′ ], no sistema de coordenadas amarrado na bicicleta.
40. Isto explica os parâmetros , , , exigidos pela função quiver(). Uma explicação esclarecedora, Mestra!
41. Voltando à expressão ′ = , . O terceiro ingrediente é a igualdade. Ela pode ser verdadeira ou não! Em outras palavras, ela expressa uma equação.
42. Não se trata de uma equação algébrica (ou transcendental) onde a incógnita é um número. Nem mesmo de um sistema de equações lineares (ou não-lineares), no qual a incógnita é um vetor. A incógnita é uma função ↦ e a expressão ′ = , relaciona e sua derivada ′ = Τ .
43. Portanto: Se uma função : ↦ () torna a igualdade ′ = (, ) verdadeira, ela é uma solução dessa equação. É nesse instante que recordo a velha dualidade: funções são vetores definidos ponto a ponto.
44. Muito oportunamente lembrada, Mestra. Assim, para que uma função ↦ () seja uma solução da equação ′ = (, ) ela precisa tornar verdadeira a igualdade ′ = , em cada ponto do seu domínio.
45. Como já vimos, um problema de valor inicial, PVI, é constituído por uma equação diferencial mais uma condição inicial: ቊ ′ = (, ) 0 = 0 . A pergunta que surge é: Todo problema de valor inicial, PVI, tem solução?
46. Pois é Loirinha, nem sempre. Alguns PVI podem possuir mais de uma solução.
47. Confira a afirmação do Mestre em “Equações Diferenciais Ordinárias”, C.L. Doering e A.O. Lopes, Teorema de Cauchy- Peano, pg.387. Para que um PVI possua solução, basta que a função (, ) ↦ , seja contínua no seu domínio e que 0 , 0 ∈ .
48. Entretanto, para a unicidade de solução de um PVI, é preciso que a função seja Lispschitz contínua no seu domínio . Para maiores detalhes veja o Teorema de Picard-Lindelöf, às pgs. 384,385 da referência anterior.
49. Livros de Matemática, muitos livros. Bons, muito bons e baratos!
50. Mestres, testar essa condição de Lipschitz é complicado. Não há algum teste mais simples? Tem sim Loirinha. Lendo, descobri no livro de Doering & Lopes que basta testarmos se e Τ (, ) são funções contínuas em (nas duas variáveis).
51. A grande pergunta é: Dada a função (, ) ↦ (, ) satisfazendo as condições de existência e unicidade, como obter, numericamente, a função ↦ (), que resolve o problema de valor inicial.
52. Nas transparências a seguir, mostrarei que existem diversas possibilidades de obter a solução numericamente.
53. Na realidade obteremos aproximações para o problema de valor inicial ቊ ′ = (, ) 0 = 0 . Começamos assumindo que a solução ↦ () é única e está definida num intervalo num intervalo [a,b]
54. Em seguida criamos uma partição , = = 0 , 1 , … , = e calculamos aproximações 0 , 1 , … , para os valores exatos 0 , 1 , … , ( ). A utilização de partições uniformes, nas quais +1 = + ℎ, para = 0,1, … , − 1 com ℎ = ( − )/, permite criar algoritmos recursivos que permitam calcular facilmente +1 a partir de
55. Conhecida a ↦ (), para ℎ pequeno, o quociente de Newton fornece ′ ≅ +1 − ℎ . Então, se ↦ () é solução da equação ′ = (, ) teremos +1 − ℎ ≅ ( , )
56. Teríamos então uma fórmula recursiva de aproximação: +1 ≅ + ℎ ∙ , , = 0,1, … − 1 Teríamos ou teremos?
57. Teríamos, porque ↦ () é nossa incógnita. Não a conhecemos! Entretanto, para = 0 a aproximação é válida, pois da condição inicial 0 = 0 tiramos 1 ≅ 0 + ℎ ∙ 0 , 0 .
58. De fato, a fórmula fornece uma aproximação 1 para (1 ): 1 ≅ 1 = 0 + ℎ ∙ 0 , 0 . Repetindo a ideia da tangente para os pontos 1 , (1 ) , 2 , (2 ) , etc, temos o método de Euler.
60. Atenção, no texto abaixo, ↦ () é a solução.
61. O Mestre fez um programa que mostra a solução do PVI ቊ ′ = (, ) 0 = 0 para ∈ [0, ] pelo método Euler, com passo ℎ escolhido pelo usuário. Além disso, comparo com a solução “exata” e mostro o erro.
62. Uma execução do programa:
63. Este é o início do programa
64. A continuação.
65. E a parte gráfica.
66. Uma outra forma de obter o método de Euler:
67. Outro método:
68. E um 3º método:
69. Este 3º método também é conhecido como método de Heun.
70. Uma interpretação geométrica, via retas tangentes, do método de Heun.
71. Este método é famoso!
72. Famoso e eficiente, é (ℎ4) no erro global.
73. Um detalhe associativo:
76. As quatro transparências a seguir são apenas informativas.
81. Nesta e nas próximas três transparências, um exemplo detalhado do método de Runge-Kutta para EDOs de 2ª ordem.
85. O que há na scipy sobre resolução de EDOs:
86. São apenas dois métodos, o odeint e o ode. Mas, conforme veremos, muito poderosos.
93. O ODE é mais que um método, trata-se de uma classe.
103. Tchau, esta foi nossa última aula!