Equacoes_diferenciais_ordinarias

Prévia do material em texto

Equações Diferenciais Ordinárias 1 
Métodos Numéricos para Equações Diferenciais Ordinárias 
 
 
Introdução 
 
 Diversos problemas técnicos e científicos são descritos matematicamente por equações 
diferenciais que representam variações das quantidades físicas que os descrevem. Alguns 
exemplos de equações diferenciais são: 
(1) Reação química de 1a ordem BA →← , descrita pela equação AA kCdt
dC
−= , na qual CA é a 
concentração do reagente A, k a constante da reação e t o tempo decorrido desde o início da 
reação. 
(2) Descarga de um circuito elétrico contendo uma resistor em série com um capacitor, descrito 
pela equação Q
C
dt
dQRV +=0 , para a qual V0 é a tensão contínua de alimentação do circuito, 
R a resistência, C a capacitância, Q a carga elétrica acumulada no capacitor e 
dt
dQi = a 
corrente do circuito. 
(3) Condução de calor num material sólido, descrito pela equação de Fourier 
dx
dTkAq =& , na qual 
q& é o fluxo térmico, k a condutividade térmica, A a área de secção transversal ao fluxo 
térmico, T a temperatura e x a coordenada espacial na direção do fluxo de calor. 
(4) Pêndulo simples, descrito pela equação θ−=θ seng
dt
d
l2
2
, na qual θ é o ângulo formado pelo 
pêndulo em relação ao eixo vertical, g a aceleração da gravidade, l o comprimento do 
pêndulo e t o tempo. 
 
 Dos exemplos citados, vemos que o grau (ou ordem) de uma equação diferencial pode 
variar. O grau de uma equação diferencial é definido pelo termo da equação que contém a 
derivada de maior ordem. Por exemplo, a seguinte equação diferencial 02 =−+xy´ é uma 
equação diferencial de 1o grau porque a derivada y´ é de 1a ordem. Já a equação diferencial 
0852 =+−+′+′′−′′′ xyyyxy é uma equação diferencial de 3o grau porque o termo de derivada 
de maior ordem é de 3a ordem. Se a solução de uma equação diferencial y for uma função de uma 
única variável x, isto é, se y = y(x), então a equação diferencial é chamada de equação diferencial 
ordinária. 
 
Definição 
 
 Uma equação diferencial ordinária de grau n é uma equação que pode ser descrita na 
forma geral como: 
 
)y,,y,y,y,x(fy )n()n( 1−′′′= K (1) 
 
sendo que 
n
n)n(
dx
ydy ≡ empregando a notação de Leibniz. 
 
Equações Diferenciais Ordinárias 2 
 Uma equação diferencial ordinária (E.D.O.) de 1a ordem para duas variáveis x e y é 
definida como uma equação da forma espacial: 
 
′ = =y
dy
dx
f x y( , ) (2) 
 
ou para duas variáveis y e t, na forma temporal como: 
 
& ( , )y dy
dt
f y t= = (3) 
 
 No caso particular, f(x,y) = f(x), podemos obter a solução geral para E.D.O. de 1a ordem 
(2) por separação de variáveis: 
 
dx)x(fdy)x(f
dx
dyy ⋅=⇒==′ (4) 
 
que pode ser integrada diretamente como: 
 
∫ +⋅= Cdx)x(fy (5) 
 
onde C é a constante de integração. Para obtermos uma solução particular (ou seja, um valor 
específico para a constante C), é necessário fornecer a condição de contorno para a equação (2): 
 
f x y C( , )0 0 0= (6) 
 
 Se y = y(x) é uma solução, então dy/dx = f(x,y) e y0 = y(x0) é a condição de contorno da 
equação (2). 
 Se considerarmos a E.D.O. (3) na qual a variável t representa o tempo, então a condição 
para obtenção de uma solução particular de (3) é chamada condição inicial (análoga à condição 
de contorno, somente que esta se aplica a problemas envolvendo apenas coordenadas espaciais). 
 
Exemplo 1 
Seja a E.D.O. de 1a ordem: yy =′ , cuja solução analítica geral é expressa por y Cex= . Se 
impusermos como condição de contorno y(0) = 1, isto é, em x = 0, y = 1 e substituirmos na 
solução geral, vem que, 1 0= =Ce C . 
Portanto, a solução particular da E.D.O. y’ = y é obtida substituindo-se o valor da constante de 
integração C calculada da condição de contorno y(0) = 1, resultando: 
 
y ex= 
 
Exemplo 2 
 Seja a E.D.O. de 1o grau, y' = x + y, cuja solução analítica, obtida pelo Método dos 
Fatores Integrantes1, é expressa por: y(x) = Cex - x - 1. Se adotarmos a condição de contorno y(0) 
 
1
 Matemática Superior, E. Kreyszig, Livros Técnicos e Científicos, Rio de Janeiro,1969, p.69. 
 
Equações Diferenciais Ordinárias 3 
= 1, vem que y(0) = C - 1 = 1. Portanto, C = 2, que substituindo na solução geral, resulta a 
solução particular: y(x) = 2ex - x - 1. 
 É importante salientar que a solução geral representa uma família de soluções (isto é, um 
conjunto infinito de soluções) e que a solução particular representa uma solução única. Como 
nos métodos numéricos pressupõe-se que a solução do problema seja única, isto irá requerer na 
descrição do problema a especificação da condição de contorno juntamente com a equação 
diferencial. 
 
 
Método de Euler 
 
 O Método de Euler é um método aproximado de 1a ordem, isto é, ele aproxima a solução 
da E.D.O. de 1o grau y(x) = y(x) por uma função de 1o grau, isto é, por uma reta. A Fig. 1 ilustra a 
aproximação da solução exata y = y(x) por uma solução aproximada y , obtida pelo 
prolongamento da reta tangente à curva de y = y(x) em x = x0 até o valor de x para o qual deseja-
se obter a solução da E.D.O. 
 
y
x
y = y(x)
Solução exata da E.D.O.
x
y0
0
Condição de
contorno
x1
y
1
y
1
Valor exato
Valor
aproximado
pelo método de
Euler
Aproximação de y(x)
pelo método de
Euler (aproximação
linear)
 
Fig. 1 Solução gráfica da E.D.O. pelo método de Euler. 
 
 A equação genérica para o cálculo da solução de uma E.D.O. de 1o grau pelo Método de 
Euler é expressa por: 
 
y y hf x yi i i i+ = +1 ( , ) (7) 
para a qual 
h x xi i= −+1 
 
Exemplo 1 
Seja a E.D.O. y’ = x + y, com a condição de contorno y(0) = 1. A solução da E.D.O. empregando 
o método de Euler será calculada no intervalo [0; 5]. 
A equação do método de Euler para a E.D.O. deste exemplo tem a forma: 
( )iiii yx.hyy ++=+1 
 
Equações Diferenciais Ordinárias 4 
(a) h = 1 
 
i = 0 x1 = x0 + h = 0 + 1 = 1 
 y1 = y0 + h.(x0 + y0) = 1 + 1.(0 + 1) = 2 
 
i = 1 x2 = x1 + h = 1 + 1 = 2 
 y2 = y1 + h.(x1 + y1) = 2 + 1.(1 + 2) = 5 
 
i = 2 x3 = x2 + h = 2 + 1 = 3 
 y3 = y2 + h.(x2 + y2) = 5 + 1.(2 + 5) = 12 
 
i = 3 x4 = x3 + h = 3 + 1 = 4 
 y4 = y3 + h.(x3 + y3) = 12 + 1.(3 + 12) = 27 
 
i = 4 x5 = x4 + h = 4 + 1 = 5 
 y5 = y4 + h.(x4 + y4) = 27 + 1.(4 + 27) = 58 
 
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0
50
100
150
200
250
300
x
y 
=
 
y(x
)
Exato
Euler
 
Fig. 2 Gráfico comparativo entre a solução exata e a solução pelo método de Euler (h = 1) 
 
 
Cálculo com h = 0,5. 
(b) h = 0,5 
Os resultados estão apresentados na tabela seguinte. 
 
i xi yi f(xi,yi) xi+1 yi+1 
0 0 1,0 1,0 0,5 1,5 
1 0,5 1,5 2,0 1,0 2,5 
2 1,0 2,5 3,5 1,5 4,25 
3 1,5 4,25 5,75 2,0 7,125 
 
Equações Diferenciais Ordinárias 5 
4 2,0 7,125 9,125 2,5 11,6875 
5 2,5 11,6875 14,1875 3,0 18,7813 
6 3,0 18,7813 21,7813 3,5 29,6719 
7 3,5 29,6719 33,1719 4,0 46,2578 
8 4,0 46,2578 50,2578 4,5 71,3867 
9 4,5 71,3867 75,8867 5,0 109,3301 
10 5,0 109,3301 
   
 
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0
50
100
150
200
250
300
x
y 
=
 
y(x
)
Exato
Euler
 
Fig. 3 Gráfico comparativo entre a solução exata e a solução pelo método de Euler (h = 0,5) 
 
Programa Matlab 
 
% euler1 Programa para o calculo da E.D.O. y' = x + y 
% Metodo de Euler 
% Condicao de contorno: y(0) = 1 
 
clear; 
% Condicao de contorno 
x(1) = 0; 
y(1) = 1; 
n = input('Numero de intervalos: '); 
xf = input('Valor de x final: '); 
h = (xf - x(1))/n; 
for i = 1:n 
 f(i) = x(i) + y(i); 
 x(i+1) = x(i) + h; 
 y(i+1) = y(i) + h*f(i); 
end 
% Calculo da solucao exataxe = 0:0.1:xf; 
ye = 2*exp(xe) - xe - 1; 
% Grafico comparativo: solucao pelo metodo de Euler e a solucao exata 
plot(xe,ye,'-r',x,y,'ob'); 
xlabel('x'); ylabel('y = y(x)'); 
legend('Exato','Euler'); shg 
 
Equações Diferenciais Ordinárias 6 
 
Exemplo 2 
Seja a E.D.O. y’ = -ky, com a condição de contorno y(1) = 1. Calcular a solução da E.D.O. 
empregando o método de Euler em x = 2, para h = 0,5 e h = 0,25. 
Neste exemplo, por questão de conveniência, vamos realizar os cálculos numa tabela que 
sumariza os resultados. 
A equação do método de Euler para a E.D.O. y’ = y é: 
 
iii yhyy .1 +=+ 
 
(a) h = 0,5 
 
i xi yi yi+1 
0 1,0 1,0 1,5 
1 1,5 1,5 2,25 
2 2,0 2,25 
 
 
(b) h = 0,25 
 
i xi yi yi+1 
0 1,0 1,0 1,25 
1 1,25 1,25 1,5625 
2 1,5 1,5625 1,9531 
3 1,75 1,9531 2,4414 
4 2,0 2,4414 
 
 
(c) A solução analítica é dada por: 
 
∫ ∫ ′+=⇒=⇒= Cxydxy
dyy
dx
dy ln 
 
Reescrevendo a solução analítica na forma y = f(x), resulta: 
 
xCey = 
 
A constante de integração C é calculada a partir da condição de contorno do problema: 
 
1111)1( −=⇒=⇒= eCCey 
 
que, substituindo na solução analítica geral, resultará na expressão: 1−= xey como solução 
analítica particular do problema. 
Calculando-se a solução exata em x = 2, obtém-se y(2) = e2-1 = e1 = 2,7183. Comparando-se o 
resultado exato com os resultados aproximados de (a) e (b), resulta: 
h = 0,5 erro = 2,7183 – 2,25 = 0,47 
h = 0,25 erro = 2,7183 – 2,4414 = 0,28 
o que corresponde a uma redução de 1,7 vezes no erro quando o intervalo h é reduzido pela 
metade. 
 
Equações Diferenciais Ordinárias 7 
Método de Euler Modificado 
 
 Como o método de Euler baseia-se no cálculo da solução pela aproximação de uma reta 
tangente traçada em x0 para avaliar a solução em x1, o Método de Euler Modificado apresenta 
uma correção na estimativa da solução em x1 calculando o valor de 1y e fazendo a média com o 
valor y0. Generalizando para o cálculo do valor estimado yi+1 , toma-se a média entre as 
inclinações das retas tangentes em xi e xi+1: 
 
( ) ( ) ( )[ ]112
1
+++= iiiiiimedio y,xfy,xfy,xf (8) 
 
na qual o argumento yi+1 = yi + h.f(xi,yi) do segundo termo de (8) é obtido a partir da solução do 
método de Euler. Substituindo o coeficiente angular médio fmedio(xi,yi) é substituído em (7): 
 
[ ]y y h f x y f x h yi i i i i i+ += + + +1 12 ( , ) ( , (9) 
 
y
x
y = y(x)
Solução exata da E.D.O.
x
y0
0
Condição de
contorno
x1
y
1
y
1
Inclinação de y(x)
em x1: correção do
método de EulerSolução em x1
pelo método de
Euler Modificado Inclinação de y(x)
em x0 pelo método
de Euler
Retas de coeficiente
angular fmedio
 
Fig. 4 Solução gráfica da E.D.O. pelo método de Euler Modificado. 
 
 
Métodos de Runge-Kutta 
 
 Os Métodos de Runge-Kutta consistem em métodos preditores-corretores de 2a e 4a 
ordem. No caso do Método de Runge-Kutta de 2a ordem, a expressão para o cálculo aproximado 
de yi+1 é equivalente à do Método de Euler Modificado, ou seja, 
 
[ ]y y h f x y f x h yi i i i i i+ += + + +1 12 ( , ) ( , (9) 
 
que pode ser re-escrito na forma: 
 
 
Equações Diferenciais Ordinárias 8 
( )
)hky,hx(fk)y,x(fk
kkhyy
iiii
ii
121
211 2
++==
++=+
 (10) 
 
 A fórmula do Método de Runge-Kutta de 4a ordem é dada por: 
 
( )
),(
)2/,2/(
)2/,2/(
),(
22
6
34
23
12
1
43211
hkyhxfk
hkyhxfk
hkyhxfk
yxfk
kkkkhyy
ii
ii
ii
ii
ii
++=
++=
++=
=
++++=+
 (11) 
 
Exemplo 
Recalculando a solução da E.D.O. y’ = x + y, com a condição de contorno y(0) = 1 no intervalo 
[0; 5] utilizando os métodos de Runge-Kutta de 2ª e 4ª ordem com h = 1, obtém-se o gráfico 
mostrado na Fig. 5. 
 
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
0
50
100
150
200
250
300
x
y 
=
 
y(x
)
Exato
RK 2a ordem
RK 4a ordem
 
Fig. 5 Gráfico comparativo entre a solução exata e as soluções pelos métodos de Runge-Kutta de 
2ª e 4ª ordem (h = 1,0). Observar que a solução aproximada pelo método de 4ª ordem é uma 
solução bastante próxima à exata. 
 
 
Programa Matlab 
 
% rungekutta Programa para o calculo da E.D.O. y' = x + y 
% pelos metodos de Runge-Kutta de 2a e 4a ordem 
% Condicao de contorno: y(0) = 1 
% Solucao analitica: y(x) = 2 exp(x) - x - 1 
 
 
Equações Diferenciais Ordinárias 9 
clear; 
x(1) = 0; % valor inicial 
y2(1) = 1; % condicao de contorno (Runge-Kutta de 2a ordem) 
y4(1) = 1; % condicao de contorno (Runge-Kutta de 4a ordem) 
h = input('Incremento h: '); 
xf = input('Valor final de x: '); 
n = floor((xf - x(1)) / h + 1); % Numero de intervalos 
for i = 1:n-1 
 x(i+1) = x(i) + h; 
 % Metodo de Runge-Kutta de 2a ordem 
 k12 = x(i) + y2(i); 
 k22 = x(i) + h + (y2(i) + h*k12); 
 y2(i+1) = y2(i) + h/2*(k12 + k22); 
 % Metodo de Runghe-Kutta de 4a ordem 
 k14 = x(i) + y4(i); 
 k24 = x(i) + h/2 + (y4(i) + h*k14/2); 
 k34 = x(i) + h/2 + (y4(i) + h*k24/2); 
 k44 = x(i) + h + (y4(i) + h*k34); 
 y4(i+1) = y4(i) + h/6*(k14 + 2*k24 + 2*k34 + k44); 
end 
% Solucao exata 
xe = 0:0.1:xf; 
ye = 2*exp(xe) - xe - 1; 
% Grafico: solucoes pelos metodos de Runge-Kutta e solucao exata 
plot(x,y2,'sr',x,y4,'ob',xe,ye,'-k'); 
xlabel('x'); ylabel('y = y(x)'); 
legend('Exato','RK 2a ordem','RK 4a ordem'); shg 
 
 
Equações Diferenciais Ordinárias de 2ª Ordem 
 
 Os métodos numéricos vistos até aqui se aplicam somente à solução de E.D.O. de 1ª 
ordem. Entretanto, uma E.D.O. de 2ª ordem pode ser escrita como um sistema de duas E.D.O. de 
1ª ordem. Assim, pode-se utilizar os métodos de Runge-Kutta na solução de E.D.O. de 2ª ordem 
e ordem superior. 
 Considerando a E.D.O. genérica de 2ª ordem, 
 
y” = f (x,y,y’) (12) 
 
que pode ser escrita como o sistema: 
 
( )z,y,xfz
zy
=′
=′
 (13) 
 
Para a qual definiu-se uma variável auxiliar z. Para o sistema de E.D.O. de 1ª ordem (13), deve-
se aplicar as condições de contorno para cada E.D.O. de 1ª ordem: y(z0) = y0 e z(x0) = z0. 
Observar que esta segunda condição de contorno se refere à condição de contorno da derivada da 
variável y = y(x). 
 
Exemplo aplicado: Pêndulo simples. 
 
 Considere o pêndulo simples mostrado na Fig. 6. 
 
 
Equações Diferenciais Ordinárias 10 
θ l
mg
m
 
 
Fig. 6 Pêndulo simples. 
 
 A equação diferencial ordinária que descreve o movimento angular θ do pêndulo é obtida 
a partir das leis de Newton: 
 
amF r
r
= (14) 
 
θθ sen2
2
mg
dt
d
m −=l (15) 
 
θθ sen
l
&& g
−= (16) 
 
 A condição inicial (condição de contorno) do problema é 0)0( θθ = (ângulo inicial do 
pêndulo). Vamos calcular agora a solução da E.D.O. de 2ª ordem (16) através do Método de 
Runge-Kutta de 2ª ordem. Como a equação (16) é uma E.D.O. de 2a ordem, precisamos convertê-
la em um sistema de E.D.O. de 1a ordem, aplicando as seguintes transformações: 
 
θ−=⇒≡θ=θ sengpp
dt
d
l
&&
 
 
de modo a obter o sistema de E.D.O. de 1a ordem: 
 
θsen
l
g
dt
dp
−= (17) 
 
p
dt
d
=
θ
 (18) 
 
com as condições iniciais: 0)0( θθ = e 0)0( pp = . 
 Para resolver este sistema de equações, calculamos a solução da equação (17) para obter o 
valor da variável p, e a equação (18) para obter a solução θ em cada intervalo de tempo t. 
 Aplicando o Método de Runge-Kutta nas equações (17) e (18), resulta, respectivamente, 
em: 
 
Equações Diferenciais Ordinárias 11 
( )
ii
ii
sen
gk,sengk
kkhpp
θ−=θ−=
++=+
ll
21
211 2
 (19) 
 
( )
ii
ii
pk,pk
kkh
==
++θ=θ +
21
211 2 (20) 
onde: 
th ∆= (21)Consideremos os seguintes valores numéricos: g = 9,8 m/s2, l = 0,5 m, θ(0) = 60º e p(0) = dθ/dt 
= 0 (velocidade inicial). 
 
0 1 2 3 4 5 6 7 8 9 10
-10
-5
0
5
10
t
p(t
) =
 
dq
/d
t (r
ad
/s
)
0 1 2 3 4 5 6 7 8 9 10
-3
-2
-1
0
1
2
3
t
q(t
) (r
ad
ia
no
s)
 
Fig. 7. Gráfico da velocidade angular p = dθ/dt e do deslocamento angular θ versus tempo 
usando o Método de Runge-Kutta de 2ª ordem (h = ∆t = 0,01 s). Observe a instabilidade da 
solução para valores crescentes de tempo. 
 
 
Equações Diferenciais Ordinárias 12 
0 1 2 3 4 5 6 7 8 9 10
-5
0
5
t
p(t
) =
 
dq
/d
t (r
ad
/s
)
0 1 2 3 4 5 6 7 8 9 10
-1.5
-1
-0.5
0
0.5
1
1.5
t
q(t
) (r
ad
ia
no
s)
 
Fig. 8. Gráfico da velocidade angular p = dθ/dt e do deslocamento angular θ versus tempo 
usando o Método de Runge-Kutta de 2ª ordem (h = ∆t = 0,001 s). Observe que esta solução é 
estável e não apresenta instabilidade na resposta para tempos crescentes. 
 
 
Programa Matlab 
 
% PENDULO Calculo da E.D.O. de 2a. ordem do pendulo simples 
% Metodo de Runge-Kutta de 2a ordem 
% Condicoes de contorno: p(0) = 0, q(0) = 60 
 
clear; 
t(1) = 0; 
p(1) = 0; 
q(1) = 60; q(1) = q(1)*pi/180; % Conversao de angulo de graus para radianos 
g = 9.81; % Aceleracao da gravidade 
L = 0.5; % Comprimento do pendulo 
h = input('Incremento h: '); 
tf = input('Valor final de t: '); 
n = floor((tf - t(1)) / h + 1); % Numero de intervalos 
for i = 1:n-1 
 t(i+1) = t(i) + h; 
 % Metodo de Runghe-Kutta de 2a ordem 
 % Calculo de p(t) 
 k11 = -g/L*sin(q(i)); 
 k21 = -g/L*sin(q(i)); 
 p(i+1) = p(i) + h/2*(k11 + k21); 
 % Calculo de q(t) 
 k12 = p(i); 
 k22 = p(i); 
 q(i+1) = q(i) + h/2*(k12 + k22); 
end 
% Graficos das solucoes de p(t) e q(t) 
figure(2); clf; 
subplot(2,1,1) 
plot(t,p); xlabel('t'); ylabel('p(t) = dq/dt (rad/s)') 
subplot(2,1,2) 
plot(t,q); xlabel('t'); ylabel('q(t) (radianos)') 
 
Equações Diferenciais Ordinárias 13 
Exercícios 
 
1. Calcular a solução das seguintes E.D.O. de 1o grau nos valores indicados, utilizando o método 
de Euler e compare com a solução exata à partir da solução analítica: 
 
(a) y’ + 2y = x2, y(0) = 0,25, y(2) 
h = 0,5 e h = 0,25 
Solução analítica: Cxxy +−=
22
2
 
 
(b) y’ + y = sen x, y(0) = -0,5, y(2) 
h = 1,0 e h = 0,5 
Solução analítica: )cos(sen xxCy −= 
 
(c) y’ + 2y = x, y(0) = 1, y(3) 
h = 1 e h = 0,5 
Solução analítica: xCexy 2
4
1
2
−+−=