Métodos numéricos para
uma EDO
resolver
Neste texto serão estudados métodos de passos simples e
passos múltiplos para resolver equação de primeira ordem.
Nos de passos simples necessitamos apenas dos resultados
de yk , do passo anterior, para determinarmos a aproximação
de yk+1. Enquanto que nos de passos múltiplos para
determinarmos a aproximação y k+1 dependemos dos valores
de y k, y k - 1 . . . .
Método de Euler
O método de Euler para resolver EDO com condições
iniciais é o método numérico mais simples. Ele consiste
em aproximar a solução y ( x ), no sentido de uma
linearização, por meio de suas tangentes, Figura 1.
Considere o problema
 dy
 f (x , y )

 dx
y ( x 0 )  y 0

ou seja, são dados um ponto de partida, (x 0 , y 0), e uma
direção a ser tomada, f ( x, y ) .
Desejamos determinar y ( z ).
A interpretação geométrica da
Figura 1 nos permite escrever a
equação:
y
y1
y(x)
y0
F ’(x 0 ) = y’ (x 0) = f (x 0 , y 0)


x
h
x0
z = x1
Interpretação geométrica do Método de Euler
Fazendo x1 – x0 = h, vamos
ter y1 = y0 + h f (x 0 , y 0) ou
F(x 1)  F(x 0) + F ’(x 0) (x1 – x0 ) (Taylor).
Diremos, portanto, que
y1  F(x 1) = F( z )
Em verdade, estamos substituindo a função desconhecida
y ( x ) por, simplesmente uma reta em todo intervalo [x 0 ; z ]
e calculando a imagem de z sobre ela o que pode ser uma
aproximação ruim para y ( z ).
Figura 1
Podemos, porém, melhorar esta
aproximação se subdividirmos o
intervalo [x 0 ; z ] em subintervalos de
amplitude constante, genericamente
chamada h, e como sabemos calcular
a direção da função incógnita y ( x )
em cada ponto, substituiremos tal
y
y(x)
y2
y1
y0
h
x0
x
h
x1
z = x2
Método de Euler considerando dois
subintervalos
Figura 2
função por um segmento de reta, em cada um destes
subintervalos. Estes segmentos terão a direção que ela
(função) tem no início de cada dos subintervalos, Figura 2.
Obtemos então:
y i + 1 = yi + h f (x i , y i ), i = 0, 1, ...
que vem a ser o método de Euler.
Exemplo : Considere o problema de valor inicial y ( 1 ) = 1 da equação
diferencial
y’ = f ( x, y ) = 2x + 3. Dividindo o intervalo [ 1; 2 ] em 1, 2 e 4 partes sucessivamente
e aplicando o método de Euler, determine o valor aproximado de y ( 2 ) para a equação
dada.
Solução:
Temos y’ = f ( x, y ) = 2x + 3, com y (1) = 1 ou seja, x 0 = 1 e y 0 = 1.
Com uma divisão do intervalo, isto é, h = 1, obtemos:
y 1 = y0 + h f (x 0 , y 0 ) = 1 + 1 [ 2 x 1 + 3 ] = 1 + 5 = 6.
Com duas divisões do intervalo, isto é, h = 0,5 , temos
y1 = y0 + h f (x 0 , y 0 ) = 1 + 0,5 [ 2 x 1 + 3 ] = 1 + 2,5 = 3,5
y2 = y1 + h f (x 1 , y 1) = 3,5 + 0,5 [ 2 x 1,5 + 3 ] = 3,5 + 3,0 = 6,5
Finalmente, considerando quatro divisões, isto é, h = 0,25, temos
y1 = y0 + h f (x 0 , y 0 ) = 1 + 0,25 [ 2 x 1 + 3 ] = 1 + 1,25 = 2,25
y2 = y1 + h f (x 1 , y 1 ) = 2,25 + 0,25 [ 2 x 1,25 + 3 ] = 2,25 + 1,375 = 3,625
y3 = y2 + h f (x 2 , y 2 ) = 3,625 + 0,25 [ 2 x 1,5 + 3 ] = 3,625+ 1,5 = 5,125
y4 = y3 + h f (x 3 , y 3 ) = 5,125 + 0,25 [ 2 x 1,75 + 3 ] = 5,125 + 1,625 = 6,75
Método Modificado de Euler
Um problema que ocorre no método “simples” de
Euler é que ele pressupõe que a função que está sendo
aproximada mantém, em todo intervalo, a direção que
ela tem no extremo “de partida” dele. O método
modificado de Euler irá considerar também uma única
direção para a função y ( x ), só que uma direção média
entre aquela do “início” do intervalo e uma estimativa
da direção no “final” dele.
Para tanto, em primeiro lugar, usando o método
“simples” de Euler, fazemos uma previsão de y i + 1,
chamada y i + 1 .
Logo,
Previsão :
y
i+1
= yi + h f (x i , y i ).
A partir desta previsão, podemos obter o valor
aproximado da direção da curva y ( x ) no ponto
( x i + 1 , y i + 1 ) através de f ( x i + 1 , y
) .
Determinamos então a chamada correção,
Correção :
i+1
__
y i + 1 = yi + h/2 [ f (x i , y i ) + f ( x i + 1 , y i + 1 )] .
Esta
expressão
é
conhecida como o método
modificado de Euler. Uma
interpretação
geométrica
deste método pode ser vista
na Figura 3.
y(x)
y
( x1 ; y1 )
Direção média
( x1 ; y1 )
x
h
x0
x1
Interpretação geométrica do
Método modificado de Euler
Figura 3
Exemplo - Encontrar a solução da
equação diferencial ordinária y’ = f ( x, y ) =
2x + 3 com a condição de valor inicial y ( 1)
= 1. Dividindo o intervalo [ 1; 2 ] em
apenas uma parte, ou seja, fazendo h =1 e,
aplicando o método de modificado de Euler,
determine o valor aproximado de y ( 2 )
para a equação dada.
Solução:
Sabendo que a cada aproximação é necessário fazer
um processo de previsão – correção e, considerando h
=1, temos
Previsão
y1
y i + 1 = yi + h f (x i , y i ) ,
no caso
= y0 + h f (x 0 , y 0 )
y = 1 + 1 f (1 ,1 ) = 1 + 1 ( 2 x 1 + 3 ) = 6
1
Correção :
y i + 1 = yi + h/2 [ f (x i , y i ) + f ( x i + 1 , y i + 1 ) ] .
y1 = 1 + 1/2 [ f (1 , 1 ) + f ( 2 , 6 ) ] =
1 + 1/2 [ 5 + 2 x 2+3 ] = 1 + 6 = 7.
Métodos de Runge-Kutta
Os métodos de Runge-Kutta são uma família de
métodos numéricos para solucionar equações
diferenciais ordinárias. São métodos que podem ser
obtidos pela série de Taylor sem a necessidade de
calcular qualquer derivada. Devido sua larga
utilização, será considerado neste texto apenas o
clássico método de Runge-Kutta de 4a ordem. Por ser
sua dedução bastante trabalhosa, limitamo-nos a
enunciar apenas sua expressão. Detalhes e provas
deste método podem ser vistas em Schwarz[25].
A expressão do método de Runge-Kutta de 4a
ordem é dada por:
y i + 1 = yi + 1/6 ( k 1 + 2 k 2 + 2 k 3 + k 4 )
,
onde
k 1 = h f ( xi , yi)
k 2 = h f ( xi + h/2 , y i + k 1 /2 )
h
k 3 = h f ( xi + /2 , y i + k 2 /2 )
k 4 = h f ( xi+ 1 , y i + k 3 )
Exemplo
- Dada a EDO a seguir,
determine o valor aproximado de y ( 1 ),
a
usando o método de Runge-Kutta de 4
ordem, considerando h=1.
dy
d x = f (x , y ) = y ;
y(0) = 1
Solução:
Usando o método de Runge-Kutta de 4a ordem, temos
1
y1 = y0 + /6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) ,
Logo,
k1= 1f(0 ,1 ) = 1
k 2 = 1 f ( 0 + 1/2 , 1 + 1/2 ) = 1 x 1,5 = 1,5
k 3 = 1 f ( 0 + 1/2 , 1 +
1,5
/2 ) = 1 x 1,75 = 1,75
k 4 = 1 f ( 1 , 1 + 1,75) = 1 x 2,75 = 2,75
1
y1 = 1 + /6 ( 1 + 2 x 1,5 + 2 x 1,75 + 2,75 )
1
= 1 + /6 ( 10,25 ) = 2,708333.
(Solução exata e = 2,7182818 )
Estendendo este problema para valores diferentes de h e
aplicando os métodos de Euler, Modificado de Euler e Runge–Kutta
a
de 4 ordem, obtemos os resultado apresentados na tabela a seguir.
No divisões
1
2
4
8
16
---
Euler
2,
2,25
2,441406
2,565784
2,637927
( 2048 div.) 2,717119
Modificado de Euler
2,5
2,640625
2,694856
2,711840
2,716590
( 128 div.) 2,718225
Runge-Kutta
2,708333
2,717346
2,718209
2,718276
2,718277
---
Métodos de passos Múltiplos
Um método de passos simples determina a aproximação yk+1
em xk+1 = xk + h usando apenas o ponto de aproximação
(xk , yk). Diferentemente deste, um método de passo múltiplo
usa as informações de valores anteriores xk-1, xk-2, . . . , xk-m
que são assumidos ser equidistantes para computar yk+1.
O Método de Adams-Bashforth / Adams-Moulton é um
método preditor - corretor como a fórmula modificada de
Euler.
A fórmula de Adams-Bashforth é dada por
yn+1 = yn + h / 24 (55yn’ - 59yn -1’ + 37yn -2’ - 9yn -3’),
onde yn’ = f(xn , yn)
Note que para calcular y4
yn -1’ = f(xn -1 , yn -1)
precisamos conhecer y0 , y1,
yn -2’ = f(xn -2 , yn -2)
y2 e y3 . O valor de y0 é a
yn-3’ = f(xn -3 , yn -3)
condição inicial e os valores
y1, y2 e y3 são calculados
por um método como a
fórmula de Runge-Kutta.
Para n  3, como preditor, podemos levar o valor de yn+1
no corretor de Adams - Moulton e assim, obtermos
yn+1 = yn + h / 24 (9yn +1’ + 19yn ’ - 5yn -1’ + yn -2’),
com yn+1’ = f (xn+1, yn+1*)
Exemplo: Aplique o método de Adams-Bashforth / AdamsMoulton com h = 0,2 para obter uma aproximação
de y(0,8) para a solução de y’ = x + y -1,
y(0) = 1.
Solução:
Com um passo h = 0,2, y(0,8) será aproximado por y4.
Inicialmente, aplicando o método de Runge-Kutta com
x0 = 0, . y0 = 1 e h = 0,2 para obtermos:
y1 = 1,0214; y2 = 1,0918
e
y3 = 1,2221.
Como f (x,y) = x + y -1, temos
y0’ = f (x0, y0) = f (0, 1) = ) + 1 - 1 = 0
y1’ = f (x1, y1) = 0,2 + 1,0214 - 1 = 0,2214
y2’ = f (x2, y2) = 0,4 + 1,0918 - 1 = 0,4918
y3’ = f (x3, y3) = 0,6 + 1,2221 - 1 = 0,8221
E assim, aplicando a fórmula de Adams-Bashforth, obtemos:
y4 = y3 + 0,2 / 24 (55y3’ - 59y2 ’ + 37y1’ - 9y0’) = 1,4254.
E para y4’ = f (x4, y4*) = 0,8 + 1,4254 -1 = 1,2254
Finalmente, y4 = y3 + 0,2 / 24 (9y4’ +19y3 ’ - 5y2’ + y1’)
= 1,4255.
(Valor exato é 1,4255)