LNCC
Métodos numéricos
para EDO’s
Prof. Paulo R. G. Bordoni
UFRJ
Slide 2
Slide 2 text
Nas transparências a seguir, mostrarei
que existem diversas possibilidades de
obter aproximações numéricas para
soluções de uma EDO.
Slide 3
Slide 3 text
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 [, ].
Slide 4
Slide 4 text
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
Slide 5
Slide 5 text
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
)
Slide 6
Slide 6 text
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
Slide 7
Slide 7 text
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
Slide 8
Slide 8 text
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
)
Slide 9
Slide 9 text
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
)
Slide 10
Slide 10 text
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).
Slide 11
Slide 11 text
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.
Slide 12
Slide 12 text
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
)
Slide 13
Slide 13 text
É. Mas o erro vai
acumulando em cada
passo.
Sim, temos o erro local ℴ(ℎ2) e
o erro acumulado − global, ℴ(ℎ).
Slide 14
Slide 14 text
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
)
Slide 15
Slide 15 text
No programa a seguir codifico o método de
Euler para aproximar o PVI
ቊ
′ = (, )
0
= 0
numa partição
, = = 0
, 1
, … ,
=
de um intervalo [, ].
Slide 16
Slide 16 text
Marquei a entrada de
dados do programa do
Mestre:
Slide 17
Slide 17 text
E aqui mostrei a
implementação do
método de Euler:
Slide 18
Slide 18 text
Esse é o gráfico da solução do
programa para o método de Euler
comparado com o resultado da
odeint( ) da Scipy:
Slide 19
Slide 19 text
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
)
ℎ
Slide 20
Slide 20 text
E nesse caso teremos a
aproximação
≅ −1
+ ℎ(
,
)
Ou, trocando por + 1:
+1
≅
+ ℎ(+1
, +1
Slide 21
Slide 21 text
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 (, ).
Slide 22
Slide 22 text
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
.
Slide 23
Slide 23 text
Vamos trabalhar uma outra
forma de obter o método de
Euler.
Slide 24
Slide 24 text
′
′()
′
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
Slide 25
Slide 25 text
Consideremos novamente uma partição
uniforme
, =
= + ℎ para = 0, ⋯ ,
Então, integrando a equação ′ = (, )
de
à +1
obtemos
+1 ′ =
+1 ,
ou
+1
−
=
+1 , .
Slide 26
Slide 26 text
A integral do lado direito pode ser
aproximada, grosseiramente, pela área azul
indicada na figura, isto é
+1 , () ≅ (
,
) × (+1
−
)
Então:
+1
−
≅ (
,
) ∙ ℎ
′
+1
′(
)
′
Slide 27
Slide 27 text
( , )
+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 !
Slide 28
Slide 28 text
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
Slide 29
Slide 29 text
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
∗ =
+ ℎ ∙ (
,
)
Slide 30
Slide 30 text
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
Slide 31
Slide 31 text
(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
Slide 32
Slide 32 text
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?
Slide 33
Slide 33 text
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
Slide 34
Slide 34 text
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( )
Slide 35
Slide 35 text
Aqui marquei o
cálculo do RK 4:
Slide 36
Slide 36 text
Aqui desenhamos
os três gráficos:
Slide 37
Slide 37 text
Os dados e a tabela
comparativa:
Slide 38
Slide 38 text
Os gráficos da solução
exata, do Runge-Kutta e
da odeint( )
Slide 39
Slide 39 text
O que há na scipy
sobre resolução de
EDOs:
Slide 40
Slide 40 text
São apenas dois métodos, o
odeint e o ode.
Mas, conforme veremos, muito
poderosos.
Slide 41
Slide 41 text
No content
Slide 42
Slide 42 text
No content
Slide 43
Slide 43 text
No content
Slide 44
Slide 44 text
No content
Slide 45
Slide 45 text
No content
Slide 46
Slide 46 text
No content
Slide 47
Slide 47 text
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.
=
Slide 48
Slide 48 text
=
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.
Slide 49
Slide 49 text
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
Slide 50
Slide 50 text
Utilizamos a . ( )
para obter a solução.
O programa está a seguir.
Slide 51
Slide 51 text
O programa.
Marquei a entrada
de dados.
Slide 52
Slide 52 text
Marquei a função que
define o sistema de eqs.
lineares e a chamada da
odeint( )
Slide 53
Slide 53 text
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.