Introdução a Resolução
Numérica de Equações
Diferenciais Ordinárias
Equações Diferenciais
Definição 6.1.1: Uma equação que envolve derivadas até
ordem n, é chamada de equação diferencial ordinária
(EDO) de ordem n e pode ser escrita na forma:
y
(n)
( x )  f ( x , y ( x ), y  ( x ),  , y
a  x b
( n  1)
( x ))
(6.1.1)
Definição 6.1.2: A solução da equação (6.1.1) é qualquer
função y = F(x) que é definida em [a,b] e tem n derivadas
neste intervalo e que satisfaz (6.1.1).
A forma mais simples de uma EDO é
y   f ( x)
(6.1.2)
onde f é contínua para a < x < b.
A solução geral desta equação é:
y( x) 

f ( x ) dx  c
(6.1.3)
com constante c determinada por y ( x 0 )  y 0
De um modo geral temos
 y n ( x )  f ( x , y , y ,  , y ( n  1 ) )


a  x  b
 y ( a )   , y  ( a )   ,  , y n 1 ( a )  
1
2
n

(6.1.4)
como por exemplo:
 y  ( x )  x y  ( x )  e x y ( x )  x 2  1  f


0  x 1
 y ( 0 )  1, y  ( 0 )  0 ,  , y  ( 0 )   1

é conveniente reduzi-la a um sistema
primeira ordem chamando : y 1  y ( x )
de EDO de
 y ( x )  x y ( x )  e x y ( x )  x 2  1  f


0  x 1
 y ( 0 )  1, y  ( 0 )  0 ,  , y  ( 0 )   1

Para reduzir a uma EDO de primeira ordem
assumimos:
y1  y ( x )
y 2  y ( x )
y 3  y ( x )
y 1 ( x )  y 2


y 2 ( x )  y 3


x
2
 y 3 ( x )  xy 2 ( x )  e y 1 ( x )  x  1  f

y 1 ( 0 )  1, y 2 ( 0 )  0 , y 3 ( 0 )   1

0 x 1
isto é, se
 y1 
~y   y 
 2
y 
 3
 y 1 


~
y    y 2 


 y 3 
 y1 (0) 


~
y ( 0)   y 2 ( 0)


 y 3 ( 0)
y2




~
~
y  
y3
  F (x, y)


 f ( x , y1 , y 2 , y 3 )
de um modo geral
~
~
 dy

 F (x, ~
y ( x ))
 dt
~
~

y
(
0
)

y0

(6.1.6)
Equações de Diferenças
Definição 6.2.1: Uma equação de diferenças de ordem n é
uma sequência de equações da forma:
 g ( y k  n , y k  n 1 , , y k )  0



yi  i

k  0 ,1,...
i  0 ,1,..., n  1
(6.2.1)
Os g k são funções de (n+1) variáveis e os valores  i são
específicos.
Definição 6.2.2: Uma solução de (6.2.1) é uma sequência
 y 0 , y 1 , y 2 , ... , y n  1 , y n , y n  1 , ... que satisfaz (6.2.1).
Uma forma especial de (6.2.1) é:
  n y k  n   n  1 y k  n  1  ...  0 y k  0

yi  i

(6.2.2)
k  0 ,1, . . .
i  0 ,1, . . . , n  1
Estas fórmulas nos lembram as fórmulas recursivas e
pode ser mostrado que toda fórmula recursiva pode ser escrita
como uma equação de diferença.
Suponha, por exemplo, que a equação (6.2.2) tenha uma
solução da forma:
yk  l
k
k=0,1,...
(6.2.5)
onde k é constante a ser determinada, então de (6.2.2) temos:
 nl
k n
  n 1 l
k  n 1
 ...  0 l  l (  n l   n  1 l
k
k
n
n 1
 ...  0 )  0
Como l  0 então p(l) = 0
p ( l )   n l   n 1 l
n
Logo se l é uma raiz de
equação (6.2.2).
n 1
 ...  0
p(l) = 0 então
(6.2.6)
satisfaz a
Exemplo :
y k  2  8 .5 y k  1  4 y k
(6.2.7)
p( l )  2 l  17 l  8  0
2
l 1  0 .5
l2  8
y k  a 1l  a 2 l
k
1
k
2
(6.2.8)
satisfaz (6.2.7). Note que, se y 0  2 e y 1  1
 a1  a 2  2

 0 .5 a 1  8 a 2  1
e a solução será

então
a1  2

a 2  0
y k  2 ( 0.5)
k
(6.2.9)
Se em vez dos valores iniciais anteriores tivermos,
por exemplo:  y 0  2 .000001

 y 1  1.000008
então teríamos
 a 1  2 .0

6
a

10
 2
, de modo que a solução
perturbada seria :
6
y k  2 ( 0.5)  10 8
k
k
(6.2.9-a)
Mesmo para valores moderados de k, a
mudança é grande.
Teorema 6.2.1: A solução da equação de diferenças
que envolve somente as raízes l 1 , l 2 , ... , l m do
polinômio característico p(l) = 0 é estável (bem
condicionada numéricamente) se todas as outras
raízes l m 1 , l m  2, ... , l n satisfazem
para k  m  1, ..., n .
l k <1
Estudaremos, agora, os métodos que aproximam uma EDO
por uma equação de diferença. Para facilitar, seja
inicialmente
 y   f ( x , y ( x ))

 y (x0 )  y0
(6.2.10)
y k 1
• Em alguns métodos para determinar
, uma
aproximação de y ( x k ) levamos em consideração apenas o
valor de f no ponto y k , isto é, temos métodos de um só passo.
• Quando estamos calculando y k  1 podemos levar em conta,
não só y k , mas também os passos anteriores y k 1 , ... , y k  l ,
etc. Tais métodos são chamados de métodos de passos
múltiplos.
A forma geral dos métodos de k passos é dada por:
k
y i 1 

j 1
k
j
y i  j 1  h   j f i  j 1
(6.2.11)
j0
Se
0  0
 Método Explícito.
Se
0  0
 Método Implícito.
Métodos de um Passo
Método de Runge Kutta (RK)
Os métodos de RK são obtidos pela série de Taylor em que se
omitem os termos de mais alta ordem na expansão. Assim,
expandindo a solução em série de Taylor obtemos:
y ( x  X )  y ( x )  X y ( x ) 
X
2
2!
y( x ) 
X
3
3!
y( x )  ...
Sabendo que
dy
df


y 



dt
dt

2
dy

f

2
 
y 
 f
 f
2


dt
x
f
x
f 

x 
x 
f
 f
t
2
f
x
 f
2
2f
xt


f
t
f f
x t
 f
2

t
2
obtemos a base para todos os métodos RK
f 
 f
y ( x  X )  y ( x )  Xf 
f




2!   x  t 
X
X 


f
3! 

3
 f
2
2
x
2
 f 
 f

 x 
2
2
(6.2.14)
2
 f 

2f


 ...
2 
xt
x t
t 

 f
2
f f
y(x)  y(xk )  y(k )


 y ( x  X )  y ( x k 1 )  y ( k  1)
 f ( y ( x ), x )  f ( y ( x ), x )  f ( k )
k
k

Discretizando, temos:
logo
f 

y ( k  1)  y ( k )  Xf ( k ) 
f (k )



2! 
x 
X
2
(6.2.17)
2
2

X  2
 f (k )

 f ( k )  

 f ( k )
   ...
 f (k )
2
 x  
3! 

x


3
levando-se em conta
f
t
 0
Método de Euler
A menor aproximação é da forma :
y ( k  1)  y ( k )  X f ( k )
k = 0,1,...
(6.2.18)
e é conhecido como MÉTODO DE EULER.
Runge Kutta de segunda ordem
Para desenvolvermos o método de segunda ordem RK2,
assumiremos que a solução assume uma aproximação
tendo a expressão:
y ( k  1)  y ( k )  c1 g 1  c 2 g 2
k = 0,1,...
(6.2.19)
onde c 1 e c 2 são constantes e
g1  Xf (k )
g 2  Xf ( y( k )   2 g1 )
com  2 constante.
(6.2.20)
(6.2.21)
Desenvolvendo em série de Taylor a expressão (6.2.21),
obtem-se:
y ( k  1)  y ( k )  ( c 1  c 2 ) X f ( k )  c 2  2 X
Logo
 c1  c 2  1

1

c  
 2 2
2

c1  c 2  1
2  1
y ( k  1)  y ( k ) 
1
2
g1  Xf (k )
2
f (k )
f ( k )
x
(6.2.23)
2 então
( g1  g 2 )
g 2  Xf ( y(k )  g1 )
 ...
Runge Kutta de quarta ordem
Normalmente X = h, por universalidade : g 1  K 1 e
g 2  K 2 , e seguindo o mesmo procedimento
obteremos RK4 (o mais importante):
y ( k  1)  y ( k ) 
1
6
( g1  2 g 2  2 g 3  g 4 )
g1  Xf (k )


 g 2  X f ( y ( k )  0 .5 g 1 )

 g 3  X f ( y ( k )  0 .5 g 2 )

 g 4  Xf ( y(k )  g 3 )
k = 0,1,...
(6.2.25)
Nos métodos RK temos as seguintes características:
 São auto-inicializáveis, isto é, a partir da
condição inicial temos condições de calcular os
demais pontos;
 Não precisam do cálculo manual de derivadas;
 Permitem fácil troca de X=h;
 São facilmente codificáveis;
 Utilizam p avaliações de função, se a ordem da
fórmula for p.
Método de Passos Múltiplos
Conforme vimos, os métodos de passo simples
precisam de informação sobre a solução apenas em
x = x n para achar uma aproximação para y( x n + h ).
A característica dos métodos de passo múltiplos é
que eles usam informações sobre a solução em mais de um
ponto. Inicialmente vamos supor que conhecemos
x 0 , x 1 , ... , x n
aproximações para y(x) em
e
, i = 0,1,... .
x i 1  x i  h
A seguir exporemos aqui uma classe de métodos de
passo múltiplo que é baseado no principio de integração
numérica (Métodos Adams-Bash Forth).
Métodos Explícitos
São obtidos quando trabalhamos com x n , x n  1 , ... , x n  m
para aproximar a integral :
xn 1
y ( x n 1 )  y ( x n ) 

f ( x , y ( x )) dx
xn
Aproximamos f ( x , y ( x )) pelo polinômio de grau m, Pm ( x )
que interpola f ( x , y ) em x n , x n  1 , ... , x n  m e então:
xn 1
y ( x n 1 )  y ( x n ) 

xn
Pm ( x ) dx
Escolhendo m = 3 então vamos usar ( x n , y n ) ; ( x n 1 , y n 1 ) ;
( x n  2 , y n  2 ) ; ( x n  3 , y n  3 ) aproximando f ( x , y ( x )) pelo
polinômio de grau 3 (  ) P3 ( x ) que interpola f ( x , y ( x ))
nos pontos acima, fazendo as contas:
y n 1  y n 
h
24
[ 55 f n  59 f n  1  37 f n  2  9 f n  3 ]
Se o polinômio fosse do quarto grau 
y n 1  y n 
h
720
[1901 f n  2774 f n  1  2616 f n  2  1274 f n  3  251 f n  4 ]
Métodos Implícitos
São obtidos quando trabalhamos com x n , x n  1 , ... , x n  m
(Adams-Moulton). Se m = 2 vamos usar
( x n 1 , y n 1 ), ( x n , y n ), ( x n 1 , y n 1 ), ( x n  2 , y n  2 )
Da mesma forma que fizemos anteriormente
x n 1
y ( x n 1 )  y ( x n )  
P3 ( x ) dx 
xn
y n 1  y n 
h
24
9 f n 1  19
f n  5 f n 1  f n  2 
Notamos para métodos de mesma ordem :
 eles usam menos informações que os métodos
explícitos;
 os coeficientes de f i são menores e portanto
são fórmulas menos sensíveis ao arredondamento;
 os limites do erro de truncamento são menores;
 o problema de cálculo de y i  1 , pois precisamos
de f i  1 , o que torna o método inaplicável
numéricamente.
Métodos Previsores Corretores
Para contornar o último quesito podemos adotar um
esquema:
1.
Calculamos o valor de y i  1 por um método
explícito e o chamamos de valor previsto.
2. Com o valor previsto, calculamos f i  1 .
3. Recalculamos y i  1 , isto é, corrigimos o valor
de y i  1 já previsto usando uma fórmula implícita.
Matlab
A biblioteca do Matlab de EDOs os seguintes métodos de
valor inicial:
ode23
método explícito de um passo, RK de ordem
baixa.
ode45
método explícito de um passo, RK de ordem
média. Este é geralmente, o primeiro
método a se tentar em um novo problema.
ode113método de passo múltiplo, de AdamsBashforth-Moulton de ordens variadas.
ODE 23
Adequado para problemas que apresentam bruscas
variações na solução, para os quais é aceitável uma baixa
precisão, ou problemas em que f(t, y) não é suave, ou seja,
descontínua.
Equação da Mola
sujeita a uma força asin (t )
w1  w 2
w 2  
c
m
w2 
h
m
w1 
a
sin ( t )
m
function wd=spring(t,w);
a=2.0; m=2.0;
c=1.4;h=0.1;
wd = [w(2);-c/m*w(2)-h/m*w(1)+a/m*sin(t)];
>> tspan=[0 100];
>> wo=[0;0];
>> [t,w] = ode23(@mola,tspan,wo);
>> plot(t,w(:,1));
>> grid
>> xlabel('tempo, segundos');
>> ylabel('Elongamento da mola, metros');
>> plot(w(:,1),w(:,2));
>> xlabel('Velocidade das oscilacoes da mola');
>> ylabel('Amplitude das oscilacoes da mola');
ODE 45
Adequado para problemas não-stiff que exijam
precisão moderada.
Equação de Van der Pol
x   (1  x ) x  x  0
2
onde  é um parâmetro positivo. Escolhendo y 1  x
e
y 2  x a equação de Van der Pol se torna:
y 1  y 2
y 2   (1  y ) y 2  y1
2
1
function yprime=vdpol(t,y)
mu=2;
yprime = [y(2);mu*(1-y(1)^2)*y(2)-y(1)];
>> tspan = [0 20];
>> yo = [2; 0];
>> [t,y] = ode45(@vdpol,tspan,yo);
>> plot(t,y(:,1),t,y(:,2),'--');
>> xlabel('tempo');
>> title('Solucao de Van der Pol');