Métodos Matemáticos para EDOs

Prévia do material em texto

Métodos Matemáticos 
Diferenças Finitas e Euler de 2ª ordem 
 
 
• Diferença finita de 2ª ordem 
 
Para soluções de EDO de 2ª ordem, precisamos definir diferenças finitas de 2ª ordem. Vamos usar 
as Séries de Taylor novamente, mas ao invés de truncar na derivada de primeira ordem, vamos 
truncar a série na derivada de 2ª ordem. 
 
𝑓(𝑥0 + ℎ)~𝑓(𝑥0) + 𝑓′(𝑥0)ℎ +
𝑓′′(𝑥0)
2
ℎ2 
 
𝑓(𝑥0 − ℎ)~𝑓(𝑥0) − 𝑓′(𝑥0)ℎ +
𝑓′′(𝑥0)
2
ℎ2 
 
Somando: 
 
𝑓(𝑥0 + ℎ) + 𝑓(𝑥0 − ℎ)~2𝑓(𝑥0) + 𝑓′′(𝑥0)ℎ2 
 
 
𝑓′′(𝑥0)~
𝑓(𝑥0 + ℎ) − 2𝑓(𝑥0) + 𝑓(𝑥0 − ℎ)
ℎ2
 
 
Código: 
 
𝑦′′~
𝑦𝑖+1 − 2𝑦𝑖 + 𝑦𝑖−1
ℎ2
 
 
Exemplo: Resolver o PVI 
 
{
𝑦′′ + 3𝑦′ + 2𝑦 = cos 𝑡
𝑦(0) = 0; 𝑦′(0) = −6 10⁄
 
 
𝑦(𝑡) = −
11
10
𝑒−𝑡 + 𝑒−2𝑡 +
1
10
cos 𝑡 +
3
10
sin 𝑡 
 
Podermos usar 3 tipos diferentes de diferenças finitas para substituir 𝑦′. Vamos resolver de duas 
formas diferentes: 
 
(i) Utilizando diferença progressiva para 𝑦′: 
 
𝑦′′ + 3𝑦′ + 2𝑦 = cos 𝑡 
 
𝑦𝑖+1 − 2𝑦𝑖 + 𝑦𝑖−1
ℎ2
+ 3 (
𝑦𝑖+1 − 𝑦𝑖
ℎ
) + 2𝑦𝑖 = cos 𝑡𝑖 
 
𝑦𝑖+1 =
ℎ2 cos 𝑡𝑖 − 𝑦𝑖(2ℎ2 − 3ℎ − 2) − 𝑦𝑖−1
1 + 3ℎ
 
 
Para a implementação do programa teremos um problema. A segunda condição inicial é dada na 
forma diferencial (𝑦′(0) = −6 10⁄ ). Para contornar essa situação, uma estratégia é usar diferença 
progressiva para estimar o valor de 𝑦(2). 
 
Suponha que as condições iniciais sejam dadas de forma genérica: 
 
{
𝑦(𝑡0) = 𝑎
𝑦′(𝑡0) = 𝑏
 
 
No programa, devemos fazer: 
 
𝑦(1) = 𝑎 
 
𝑏 = 𝑦′ =
𝑦(2) − 𝑦(1)
ℎ
 ⇒ 𝑦(2) = 𝑦(1) + 𝑏ℎ 
 
Matlab 
 
f = @(t) -11*exp(-t)/10 + exp(-2*t) + cos(t)/10 + 3*sin(t)/10; 
 
a = 0; 
b = -6/10; 
 
y(1) = a; 
t(1) = 0; 
 
h = 0.1; 
 
y(2) = y(1) + b*h; % Dif. Prog. 
t(2) = t(1) + h; 
 
for i=2:80 % 8 segundos 
 y(i + 1) = (h^2*cos(t(i)) - y(i)*(2*h^2 - 3*h - 2) - y(i - 1))/(1 + 3*h); 
 t(i + 1) = t(i) + h; 
end 
 
plot(t,y,'*r') 
hold on 
plot(t,f(t),'b') 
title('Dif. Finitas - 2a Ordem - PROGRESSIVA') 
 
 
Note que o resultado inicial não é muito bom. Podemos ter um resultado alternativo utilizando 
diferença centrada ao invés da progressiva. 
 
(ii) Utilizando diferença centrada para 𝑦′: 
 
𝑦′′ + 3𝑦′ + 2𝑦 = cos 𝑡 
 
𝑦𝑖+1 − 2𝑦𝑖 + 𝑦𝑖−1
ℎ2
+ 3 (
𝑦𝑖+1 − 𝑦𝑖−1
2ℎ
) + 2𝑦𝑖 = cos 𝑡𝑖 
 
𝑦𝑖+1 =
ℎ2 cos 𝑡𝑖 − 𝑦𝑖(2ℎ2 − 2) − 𝑦𝑖−1(1 − 1,5ℎ)
1 + 1,5ℎ
 
 
Matlab 
 
f = @(t) -11*exp(-t)/10 + exp(-2*t) + cos(t)/10 + 3*sin(t)/10; 
 
a = 0; 
b = -6/10; 
 
y(1) = a; 
t(1) = 0; 
 
h = 0.1; 
 
y(2) = y(1) + b*h; % Dif. Prog. 
t(2) = t(1) + h; 
 
for i=2:80 % 8 segundos 
 y(i + 1) = (h^2*cos(t(i)) - y(i)*(2*h^2 - 2) - y(i - 1)*(1 - 1.5*h))/(1 + 1.5*h); 
 t(i + 1) = t(i) + h; 
end 
 
plot(t,y,'*r') 
hold on 
plot(t,f(t),'b') 
title('Dif. Finitas - 2a Ordem - CENTRADA') 
 
 
 
• Euler de 2ª ordem 
 
 
Exemplo: Resolver o PVI 
 
{
𝑦′′ + 3𝑦′ + 2𝑦 = cos 𝑡
𝑦(0) = 0; 𝑦′(0) = −6 10⁄
 
 
𝑦(𝑡) = −
11
10
𝑒−𝑡 + 𝑒−2𝑡 +
1
10
cos 𝑡 +
3
10
sin 𝑡 
 
Estratégia: 2 EDOs de 1ª ordem nas variáveis 𝑦 (original) e 𝑧 (variável nova) 
 
Mudança de variável: 
 
{
𝑦′ = 𝑧
𝑧′ = 𝑦′′ 
 
EDO Função 𝒇(𝒕, 𝒚, 𝒛) Código 
𝑦′ = 𝑧 𝑓1(𝑡, 𝑦, 𝑧) = 𝑧 𝑦𝑛+1 = 𝑦𝑛 + 𝑓1(𝑡, 𝑦, 𝑧)ℎ 
𝑧′ = 𝑦′′ = cos 𝑡 − 2𝑦 − 3𝑧 𝑓2(𝑡, 𝑦, 𝑧) = cos 𝑡 − 2𝑦 − 3𝑧 𝑧𝑛+1 = 𝑧𝑛 + 𝑓2(𝑡, 𝑦, 𝑧)ℎ 
 
 
Matlab 
 
% Condições iniciais nas variáveis y e z 
y(1) = 0; 
z(1) = -6/10; % -6/10 = y'(1) = z(1) 
 
t(1) = 0; 
h = 0.1; 
 
f1 = @(t,y,z) z; 
f2 = @(t,y,z) cos(t) - 2*y - 3*z; 
 
for i=1:50 
 y(i + 1) = y(i) + f1(t(i),y(i),z(i))*h; 
 z(i + 1) = z(i) + f2(t(i),y(i),z(i))*h; 
 t(i + 1) = t(i) + h; 
end 
 
plot(t,y,'*r') 
hold on 
 
ya = @(t) -11*exp(-t)/10 + exp(-2*t) + cos(t)/10 + 3*sin(t)/10; 
plot(t,ya(t),'b') 
 
 
Exercício (1,0 ponto na P3. Entrega na data da P3): Considere o PVI 
 
{
𝑦′′ + 2𝑦′ + 𝑦 = sin 𝑡
𝑦(0) = 2; 𝑦′(0) = −1
 
 
𝑦(𝑡) =
(7𝑡 + 5)𝑒−𝑡 − cos 𝑡
2
 
 
(a) Plote em uma mesma janela a solução analítica e a solução numérica pelo método de Euler 
 
(b) Plote em uma nova janela, a solução analítica e a solução numérica pelo método das diferenças 
finitas, utilizando a diferença centrada.