Aula 6 - EDO

Prévia do material em texto

Computação Numérica
Aula de Laboratório 6: 
Solução de Equações Diferenciais Ordinárias (EDOs)
Professor Filipe Taveiros
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑𝑦
𝑑𝑥
= 𝑓(𝑥)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑𝑦
𝑑𝑥
= 𝑓(𝑥, 𝑦)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑2𝑦
𝑑𝑥2
= 𝑓(𝑥)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑2𝑦
𝑑𝑥2
= 𝑓(𝑥, 𝑦)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑2𝑦
𝑑𝑥2
= 𝑓(𝑥, 𝑦, 𝑦′)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑2𝑦
𝑑𝑥2
= 𝑓(𝑥, 𝑦, 𝑦′)
Aquecimento
Derivada primeira
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑𝑛𝑦
𝑑𝑥𝑛
= 𝑓(𝑥)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑𝑛𝑦
𝑑𝑥𝑛
= 𝑓(𝑥, 𝑦)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑𝑛𝑦
𝑑𝑥𝑛
= 𝑓(𝑥, 𝑦, 𝑦′)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑𝑛𝑦
𝑑𝑥𝑛
= 𝑓(𝑥, 𝑦, 𝑦′, 𝑦′′)
Aquecimento
O problema da EDO é basicamente o mesmo da
integração que vimos na aula passada, com a
diferença de que a função a ser integrada é
expressa numa forma padrão (ordinária) em
função das suas próprias derivadas (diferencial).
𝑑𝑛𝑦
𝑑𝑥𝑛
= 𝑓(𝑥, 𝑦, 𝑦′, 𝑦′′, 𝑦′′′, … , 𝑦𝑛−1)
Aquecimento
Outra diferença é que com EDOs geralmente
estamos interessados na forma de onda da função
integrada (anti-derivada ou primitiva), e não na
integral definida (área sob a curva).
Exemplo:
Aquecimento
𝑖
Outra diferença é que com EDOs geralmente
estamos interessados na forma de onda da função
integrada (anti-derivada ou primitiva), e não na
integral definida (área sob a curva).
Exemplo:
Aquecimento
𝑖
Para ajudar a imaginação:
Em quanto tempo o 
celular vai carregar??
Outra diferença é que com EDOs geralmente
estamos interessados na forma de onda da função
integrada (anti-derivada ou primitiva), e não na
integral definida (área sob a curva).
Exemplo:
Aquecimento
−𝑉𝑖𝑛 + 𝑉𝑅 + 𝑉𝑐 = 0
𝑖
Outra diferença é que com EDOs geralmente
estamos interessados na forma de onda da função
integrada (anti-derivada ou primitiva), e não na
integral definida (área sob a curva).
Exemplo:
Aquecimento
−𝑉𝑖𝑛 + 𝑅𝑖 + 𝑉𝑐 = 0
𝑖
Outra diferença é que com EDOs geralmente
estamos interessados na forma de onda da função
integrada (anti-derivada ou primitiva), e não na
integral definida (área sob a curva).
Exemplo:
Aquecimento
𝑖
−𝑉𝑖𝑛 + 𝑅𝐶
𝑑𝑉𝑐
𝑑𝑡
+ 𝑉𝑐 = 0
Outra diferença é que com EDOs geralmente
estamos interessados na forma de onda da função
integrada (anti-derivada ou primitiva), e não na
integral definida (área sob a curva).
Exemplo:
Aquecimento
𝑖
𝑑𝑉𝑐
𝑑𝑡
=
1
𝑅𝐶
(𝑉𝑖𝑛 − 𝑉𝑐)
EDO:
𝑑𝑉𝑐
𝑑𝑡
= 𝑓 𝑡, 𝑉𝑐 =
1
𝑅𝐶
(𝑉𝑖𝑛 − 𝑉𝑐)
Condições iniciais: 𝑉𝐶 𝑡 = 0 = 𝑣0 = 0
Parâmetros 𝑅 = 𝐶 = 1 𝑉𝑖𝑛 = 12
Solução exata: 𝑉𝑐 𝑡 = 𝑒
−
𝑡
𝑅𝐶 (𝑣0 − 𝑉𝑖𝑛) + 𝑉𝑖𝑛
Aquecimento
Vamos aquecer!!
1. Defina a função 𝑓 𝑡, 𝑉𝑐 no Scilab
2. Também defina a solução analítica de
𝑉𝑐 𝑡 , vamos usá-la como referência.
3. Compare a solução analítica com os
métodos numéricos estudados.
Aquecimento
Aquecimento
Scilab 5.5.1 Console
𝑓 𝑡, 𝑉𝑐 =
1
𝑅𝐶
(𝑉𝑖𝑛 − 𝑉𝑐) 𝑉𝑐 𝑡 = 𝑒
−
𝑡
𝑅𝐶 (𝑣0 − 𝑉𝑖𝑛) + 𝑉𝑖𝑛
Aquecimento
-->R=1; C=1; vin=12; v0=0;
Scilab 5.5.1 Console
𝑓 𝑡, 𝑉𝑐 =
1
𝑅𝐶
(𝑉𝑖𝑛 − 𝑉𝑐) 𝑉𝑐 𝑡 = 𝑒
−
𝑡
𝑅𝐶 (𝑣0 − 𝑉𝑖𝑛) + 𝑉𝑖𝑛
Aquecimento
-->R=1; C=1; vin=12; v0=0;
-->deff('dvc=f(t,vc)','dvc=(1/R*C)*(vin-vc)');
Scilab 5.5.1 Console
𝑓 𝑡, 𝑉𝑐 =
1
𝑅𝐶
(𝑉𝑖𝑛 − 𝑉𝑐) 𝑉𝑐 𝑡 = 𝑒
−
𝑡
𝑅𝐶 (𝑣0 − 𝑉𝑖𝑛) + 𝑉𝑖𝑛
Aquecimento
-->R=1; C=1; vin=12; v0=0;
-->deff('dvc=f(t,vc)','dvc=(1/R*C)*(vin-vc)');
-->deff('y=vc(t)', 'y=exp(-t/(R*C))*(v0-vin)+vin')
Scilab 5.5.1 Console
𝑓 𝑡, 𝑉𝑐 =
1
𝑅𝐶
(𝑉𝑖𝑛 − 𝑉𝑐) 𝑉𝑐 𝑡 = 𝑒
−
𝑡
𝑅𝐶 (𝑣0 − 𝑉𝑖𝑛) + 𝑉𝑖𝑛
Métodos estudados
Métodos estudados
1. Euler (RK1)
Métodos estudados
1. Euler (RK1)
2. Runge-Kutta 2ª Ordem (RK2)
Métodos estudados
1. Euler (RK1)
2. Runge-Kutta 2ª Ordem (RK2)
3. Runge-Kutta 4ª Ordem (RK4)
Euler
hyxfyy
hxx
iiii
ii




),(1
1
valor da EDO no ponto
(xi ,yi)
Intervalo entre os pontos (xi ,yi) e (xi+1 ,yi+1) 𝒉 =
𝒃−𝒂
𝑵
Euler
function [x,y] = euler(f,a,b,y0,h)
x=a:h:b;
y(:,1)=y0;
for i=2:length(x)
y(:,i) = y(:,i-1) + h*f(x(i-1), y(:,i-1));
end
endfunction
Scinotes
hyxfyy
hxx
iiii
ii




),(1
1
Dessa vez passaremos a função 
como argumento de entrada
É mais comum usar o intervalo 
entre pontos h ao invés do número
de subintervalos N
Euler
--> [x,y]=euler(f,0,10,0,1);
--> plot(x,y)
Scilab 5.5.1 Console
Euler
--> [x,y]=euler(f,0,10,0,.5);
--> plot(x,y)
Scilab 5.5.1 Console
Euler
--> [x,y]=euler(f,0,10,0,.25);
--> plot(x,y)
Scilab 5.5.1 Console
Euler
--> [x,y]=euler(f,0,10,0,.125);
--> plot(x,y)
Scilab 5.5.1 Console
RK2
 
2
211
h
kkyy ii 

k1  f (x i,yi)
k2  f (xi  h,yi  k1h)
Com:
RK2
function [x,y] = rk2(f,a,b,y0,h)
x=a:h:b;
y(:,1)=y0;
for i=2:length(x)
k1 = f(x(i-1), y(:,i-1));
k2 = f(x(i), y(:,i-1)+h*k1);
y(:,i) = y(:,i-1) + (h/2)*(k1+k2);
end
endfunction
Scinotes
 
2
211
h
kkyy ii 

k1  f (x i,yi)
k2  f (xi  h,yi  k1h)
RK2
--> [x,y]=rk2(f,0,10,0,1);
--> plot(x,y)
Scilab5.5.1 Console
RK2
--> [x,y]=rk2(f,0,10,0,.5);
--> plot(x,y)
Scilab 5.5.1 Console
RK2
--> [x,y]=rk2(f,0,10,0,.25);
--> plot(x,y)
Scilab 5.5.1 Console
RK4
),(
2
1
,
2
1
2
1
,
2
1
),(
3412
231
hkyhxfkhkyhxfk
hkyhxfkyxfk
iiii
iiii















yi1  yi 
1
6
(k1  2k2  2k3  k4 )h
RK4
function [x,y] = rk4(f,a,b,y0,h)
x=a:h:b; y(:,1)=y0;
for i=2:length(x)
k1 = f(x(i-1), y(:,i-1));
k2 = f(x(i-1)+h/2, y(:,i-1)+h/2*k1);
k3 = f(x(i-1)+h/2, y(:,i-1)+h/2*k2);
k4 = f(x(i-1)+h, y(:,i-1)+h*k3);
y(:,i) = y(:,i-1) + (h/6)*(k1+2*k2+2*k3+k4);
end
endfunction
Scinotes

yi1  yi 
1
6
(k1  2k2  2k3  k4 )h
),(
2
1
,
2
1
2
1
,
2
1
),(
3412
231
hkyhxfkhkyhxfk
hkyhxfkyxfk
iiii
iiii














RK4
--> [x,y]=rk4(f,0,10,0,1);
--> plot(x,y)
Scilab 5.5.1 Console
RK4
--> [x,y]=rk4(f,0,10,0,.5);
--> plot(x,y)
Scilab 5.5.1 Console
RK4
--> [x,y]=rk4(f,0,10,0,.25);
--> plot(x,y)
Scilab 5.5.1 Console
E se a EDO mudar em algum 
momento?
• Considere que retiramos a bateria do circuito 
no instante 𝑡 = 10. A EDO torna-se:
𝑑𝑉𝑐
𝑑𝑡
= 𝑓 𝑡, 𝑉𝑐 =
1
𝑅𝐶
(−𝑉𝑐)
-->deff('dvc=f(t,vc)','if t<10 dvc=(1/R*C)*(vin-vc); else
dvc=(1/R*C)*(-vc); end');
Scilab 5.5.1 Console
RK4
--> [x,y]=rk4(f,0,20,0,.25);
--> plot(x,y)
Scilab 5.5.1 Console
Aplicação para sistemas de 
ordem superior
• Se colocássemos um indutor no circuito?
𝐿𝐶
𝑑2𝑉𝑐
𝑑𝑡2
+ 𝑅𝐶
𝑑𝑉𝑐
𝑑𝑡
+ 𝑉𝑐 = 𝑉𝑖𝑛
• Não podemos aplicar os métodos
a sistemas de ordem N, mas podemos
aplicar a N sistemas de ordem 1.
Aplicação para sistemas de 
ordem superior
• Fazendo uma mudança de variáveis
𝑉𝑐 = 𝑧1
𝑉𝑐
′ = 𝑧2
• E olhando para 
𝐿𝐶
𝑑2𝑉𝑐
𝑑𝑡2
+ 𝑅𝐶
𝑑𝑉𝑐
𝑑𝑡
+ 𝑉𝑐 = 𝑉𝑖𝑛
• Podemos substituir:
𝑧1
′ = 𝑧2
𝑧2
′ =
1
𝐿𝐶
(𝑉𝑖𝑛−𝑧1 − 𝑅𝐶𝑧2)
Aplicação para sistemas de 
ordem superior
• Assim:
𝑧 =
𝑧1
𝑧2
𝑑𝑧
𝑑𝑡
=
𝑧1
′
𝑧2
′ =
𝑧2
1
𝐿𝐶
(𝑉𝑖𝑛−𝑧1 − 𝑅𝐶𝑧2)
• Portanto
𝑓 𝑡, 𝑧 =
𝑧2
1
𝐿𝐶
(𝑉𝑖𝑛−𝑧1 − 𝑅𝐶𝑧2)
𝑧 0 =
𝑧1(0)
𝑧2(0)
Aplicação para sistemas de 
ordem superior
• Assim:
𝑧 =
𝑧1
𝑧2
𝑑𝑧
𝑑𝑡
=
𝑧1
′
𝑧2
′ =
𝑧2
1
𝐿𝐶
(𝑉𝑖𝑛−𝑧1 − 𝑅𝐶𝑧2)
• Portanto
𝑓 𝑡, 𝑧 =
𝑧2
1
𝐿𝐶
(𝑉𝑖𝑛−𝑧1 − 𝑅𝐶𝑧2)
𝑧 0 =
𝑧1(0)
𝑧2(0)
Vetor de funções
Vetor de condições iniciais
--> L=1; //Considerando L = 1
--> deff('y=f(t,z)','y=[z(2); (1/(L*C))*(12-z(1)-R*C*z(2))]')
--> [x,z]=rk4(f,0,10,[0;0],.25);
--> plot(x,z(1,:))
Scilab 5.5.1 Console
Aplicação para sistemas de 
ordem superior
𝑓 𝑡, 𝑧 =
𝑧2
1
𝐿𝐶
(𝑉𝑖𝑛−𝑧1 − 𝑅𝐶𝑧2)
--> L=1; //Considerando L = 1
--> deff('y=f(t,z)','y=[z(2); (1/(L*C))*(12-z(1)-R*C*z(2))]')
--> [x,z]=rk4(f,0,10,[0;0],.25);
--> plot(x,z(1,:))
Scilab 5.5.1 Console
Aplicação para sistemas de 
ordem superior
𝑓 𝑡, 𝑧 =
𝑧2
1
𝐿𝐶
(𝑉𝑖𝑛−𝑧1 − 𝑅𝐶𝑧2)
Vetor de funções
Vetor de condições iniciais
--> L=1; //Considerando L = 1
--> deff('y=f(t,z)','y=[z(2); (1/(L*C))*(12-z(1)-R*C*z(2))]')
--> [x,z]=rk4(f,0,10,[0;0],.25);
--> plot(x,z(1,:))
Scilab 5.5.1 Console
Aplicação para sistemas de 
ordem superior
Em qual dos z’s estamos interessados?
𝑉𝑐 = 𝑧1
𝑉𝑐
′ = 𝑧2
𝑓 𝑡, 𝑧 =
𝑧2
1
𝐿𝐶
(𝑉𝑖𝑛−𝑧1 − 𝑅𝐶𝑧2)
Vetor de funções
Vetor de condições iniciais
Aplicação para sistemas de 
ordem superior