TRABALHO EDOs e EDPsL Elizete Rocha

Prévia do material em texto

Trabalho Final: EDP’s.
Elizete Rocha da Costa
01 de julho de 2015
Relato´rio Trabalho Final: EDP’s.
1 Objetivo
Este relato´rio tem por objetivo mostrar a resoluc¸a˜o de alguns problemas de valor ini-
cial (P.V.I) e problemas de valor de contorno (P.V.C) por meio de diferentes me´todos
nume´ricos expl´ıcitos e impl´ıcito. A equac¸a˜o diferencial y˙ = −y e´ utilizada para compa-
rar os me´todos de Euler, e Runge-Kutta de quarta ordem. Para resolver o problema de
valor de contorno utilizou-se os metodos de shooting e equil´ıbrio. O sistema de equac¸o˜es
diferenciais ordina´rias referente ao exerc´ıcio 4, cuja equac¸a˜o y˙(t) + a(t) = 0, y(0) = y0
Onde a constante a e´ uma constante positiva e e´ resolvido com o me´todo de Runge-
Kutta de ordem 4. O sistema de equac¸o˜es diferenciais ordina´rias referentes ao problema
de pendulo simples e´ resolvida, utilizando o me´todo de Euler explicito e Euler impl´ıcito
2 Introduc¸a˜o
Fenoˆmenos naturais comumente sa˜o modelados na engenharia, e sa˜o descritos pelas
Equac¸o˜es diferenciais ordina´rias (EDO) , e comum se deparar com situac¸o˜es em que
muitos dos problemas, sa˜o modelados com equac¸o˜es de problema de valor inicial (PVI)
ou de valor contornos (PVC) , as quais na˜o possuem soluc¸a˜o anal´ıtica ou sa˜o muito
complexos. A utilizac¸a˜o de me´todos nume´ricos na resoluc¸a˜o de equac¸o˜es diferenciais
permite a obtenc¸a˜o de soluc¸o˜es aproximadas dentro da faixa de interesse, tornando as-
sim a resoluc¸a˜o facilitada. Neste relato´rio sera´ apresentado os procedimentos nume´ricos
para resoluc¸a˜o de EDO’s. Onde sera´ abordado os me´todos de Euler e Runge-Kutta, em
alguns exerc´ıcios, onde se fara´ uma pequena comparac¸a˜o de erro. O problema de valor
inicial, de primeira ordem pode ser escrito a partir da equac¸a˜o diferencial 1.
dy
dt
= f (t, y) (1)
y (tn+1) = y (tn) +
∫ tn+1
tn
f (t, y (t)) dt (2)
E da condic¸a˜o inicial;
y (t0) = y0 (3)
O me´todo de Euler e´ o me´todo mais simples de resolver problemas de EDO’s. Este
me´todo e´ expresso pela equac¸a˜o:
1
Figura 1: Me´todo de Euler.
yn+1 = yn + hnf (tn, yn) (4)
onde h denota o intervalo (tn+1 − tn). O me´todo de Euler consiste em calcular repetida-
mente a soluc¸a˜o do pro´ximo ponto utilizando o valor da taxa de variac¸a˜o calculada no
ponto anterior. Para se implementar o algoritmo de Euler, cuja a estrutura e´ mostrada
abaixo, a qual pode ser implementada em uma linguagem de programac¸a˜o.
Defina a funcao dy/dt = f(t,y)
Defina a condicao inicial
t(1) = t_0
y(1) = y_0
Para n=1:n-1
h = t(n+1)-t(n)
y(n+1) = y(n) + h*f(t(n),y(n))
Fim
Uma melhor aproximac¸a˜o pode ser obtida usando o me´todo de Euler impl´ıcito ou
equacao Heun, tomando a me´dia de f em ambos os pontos extremos:
yn+1 = yn +
f (tk, yn) + f (tn+1, yn+1)
2
(5)
Como o termo yn + 1 no lado direito da equac¸a˜o 5 na˜o e´ conhecido a priore, yn+1
e´ definido implicitamente. Uma expressa˜o explicita pode ser obtida aproximando yn+1
pelo me´todo de Euler:
yn+1 = yn +
f (tn, yn) + f (tn+1, yn + f (tn, yk))
2
(6)
A equac¸a˜o 6 e´ conhecida como me´todo de Heun ou me´todo de Runge-Kutta de
segunda ordem.
Tais me´todos como, de Euler, Heun pertecem tambe´m e´ conhecido como me´todos
de Runge-Kutta de primeira e segunda ordem. O mais utilizado e´ o me´todo de quarta
ordem que e´ o mais conhecido. O me´todo de Runge-Kutta de ordem 4 e´ dado por:
yn+1 = yn +
h
6
(K1 + 2 ∗K2 + 2 ∗K3 +K4) (7)
onde,
K1 = f (tn, yn) (8)
K2 = f
(
tn +
h
2
, yn +
h
2
K1
)
(9)
K3 = f
(
tn +
h
2
, yn +
h
2
K2
)
(10)
K4 = f (tn + h, yn + h ∗K3) (11)
2
Os demais n termos da equac¸a˜o 7 pode ser interpretado como uma inclinac¸a˜o me´dia.
Onde, K1 representando a derivada em tn, K2 e´ a derivada no ponto me´dio utilizando o
me´todo de Euler, K3 fornece uma segunda estimativa para a declividade no ponto me´dio
eK4 e´ o valor da inclinac¸a˜o no ponto tn+1 utilizando o me´todo de Euler. Tais me´todos sa˜o
classificados como me´todos a passos mu´ltiplos, cuja aproximac¸a˜o da resoluc¸a˜o anal´ıtica
dependera´ da quantidade de passos utilizados.
3 EXERCI´CIOS
3.1 Exerc´ıcio no 1
Implemente os algoritmos de Euler expl´ıcito e impl´ıcito para o problema do peˆndulo
simples 2 linearizado.
Figura 2: Pendulo Simples Ex.no 1
Neste exerc´ıcio sera´ mostrado a resoluc¸a˜o anal´ıtica, a resoluc¸a˜o com os me´todos de
Euler expl´ıcito e impl´ıcito para a resoluc¸a˜o do problema do peˆndulo. Para a modelagem
matema´tica do movimento de um peˆndulo simples de massa m e comprimento l e´ descrito
pela func¸a˜o θ(t) que satisfaz a equac¸a˜o diferencial 12
dθ2
dt2
+
g
l
sen (θ) = 0, (12)
f (θ) = −g
v
cos (θ) (13)
Onde o θ(t) e´ a func¸a˜o inco´gnita, onde tetha e´ a varia´vel dependente e t e´ a varia´vel
independente.
3.2 Exerc´ıcio no 2
Aplique o me´todo de Euler para o problema de valor inicial mostrados em 14 e 15 , com
h = 0.25. Compare o resultado com a soluc¸a˜o exata y(t) = e−t Repita os ca´lculos para
h/2 e para h/4. Resolva o mesmo P.V.I. utilizando o me´todo de Runge-Kutta de quarta
ordem. Comente os resultados.
Neste exerc´ıcio inicialmente, foi realizado a implementac¸a˜o de tres me´todos e reali-
zados as ana´lises Os me´todos de Euler, Heun, de 2a ordem e Runge-Kutta de 4a ordem,
que foram utilizados na implementac¸a˜o dos algoritmos para resoluc¸a˜o do problema de
valor inicial dado pela equac¸a˜o diferencial 14 com valor inicial em 15:
dy
dt
= −y, 0 ≤ t ≤ 1 (14)
t (0) = 0, y (0) = 1; (15)
As todas rotinas realizadas na resoluc¸a˜o do exerc´ıcio 2, foram implementas em Ma-
tlab, sa˜o mostradas nos anexos 1. Os resultados obtidos com a implementac¸a˜o dos
co´digos de Euler, com espac¸amentos h = 0.25, h = 0.125 e h = 0.0625 sa˜o mostrados nas
figuras 3, 4 e 5 respectivamente. Nas figuras tambe´m pode ser observado a comparac¸a˜o
3
dos valores das resoluc¸o˜es anal´ıticas, para cada valor do passo. Pode se observar que, a
medida que o passo e´ diminu´ıdo, o erro diminui consideravelmente.
Pode se observar que o me´todo de Euler, apesar de ser de fa´cil implementac¸a˜o, para
se ter uma aproximac¸a˜o do valor da soluc¸a˜o anal´ıtica, e´ necessa´rio, implementar co´digos
com passos bem pequenos. Nos resultados de cada co´digo pode se verificar que o erro
diminui a medida que a largura do intervalo diminui. Para um passo de 0,025, apresentou,
um Erro da Integrac¸a˜o Nume´rica de Euler de 19.46 por cento. Quando o h e´ diminu´ıdo
para 0,125, o erro da Integrac¸a˜o Nume´rica de Euler ja´ ca´i para 7.60 por cento. E quando
esse numero de intervalo aumenta para 16, tendo um h = 0,0625. O Erro da Integrac¸a˜o
Nume´rica de Euler ja´ e´ de 3.43 por cento. Apesar de simples a implementac¸a˜o ha´ uma
quantidade maior de ca´lculos para a aproximac¸a˜o da resoluc¸a˜o anal´ıtica. Tais observac¸o˜es
podem ser comprovadas nas figuras 3, 4 e 5.
Fazendo a mesma resoluc¸a˜o so´ que utilizando, 7. passa-se a utilizar me´todos de or-
dem superiores como o me´todo de Runge-Kutta de 4a ordem, obtems-e uma aproximac¸a˜o
melhor da soluc¸a˜o anal´ıtica, usando o mesmo conjunto de intervalos. A seguir e´ mos-
trado a os gra´ficos gerados pela implementac¸a˜o dos co´digos, usando me´todos de passos
mu´ltiplos, onde pode-se fazer a comparac¸a˜o entre soluc¸a˜o anal´ıtica, me´todo de Euler, e
de Runge-Kutta de 4a ordem.
As figuras 6, 7 e 8 mostram a aproximac¸a˜o com menor numero de passos implemen-
tando algoritmo usando me´todos de ordens elevada. Quando se tem poucos intervalos o
me´todo cla´ssico, se mostra mais eficiente. Pore´m quando a quantidade de intervalos sa˜o
maiores, essas discrepaˆncia e´ diminu´ıda. Todavia me´todo de Euler, de primeira ordem,
e´ o que mais se distancia da soluc¸a˜o anal´ıtica. Dentre os me´todos de segunda ordem,
o me´todo de Heun apresenta uma melhor concordaˆnciaem relac¸a˜o a` soluc¸a˜o exata. O
me´todo cla´ssico de Runge-Kutta de quarta ordem apresenta a melhor concordaˆncia den-
tre os me´todos testados e sua curva se sobrepo˜e a` soluc¸a˜o anal´ıtica para os quatro valores
de h analisados, como pode-se conferir na figura6.
4
3.3 Exerc´ıcio 3
4 Anexo 1 - Rotinas Exerc´ıcios 2
4.1 Rotina me´todo de Euler para h = 0,25
% Exercicio Lista 5 ELIZETE ROCHA
% Me´todo de Euler para integrac¸~ao nume´rica
% Exercicio 2 y’ = -y => y = e^x para h = 0,25
clear all
close all
clc
x0 = 0;
xf = 1;
y0 = 1;
n = 4;
deltaX = (xf - x0)/(n-1);
%y’ = y => y = e^x
x(1) = x0;
y(1) = y0;
for i = 1: n-1
y(i+1) = y(i) + deltaX * -y(i);
x(i+1) = x(i) + deltaX;
end
y_Ex =exp(-x);
er = (y(end) - y_Ex(end)) / y_Ex(end);
Erro =er* (-100);
disp(sprintf(’============================================================’))
disp(sprintf(’Valor exato da Integral e´: %0.5f ’, y_Ex(end)))
disp(sprintf(’============================================================’))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Euler: %0.5f ’, y(end)))
disp(sprintf(’============================================================’))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Euler e´: %0.2f por centoErro))
disp(sprintf(’============================================================’))
plot(x,y,x,y_Ex,’r--’,’LineWidth’,2)
xlim([0,1]);
title(’\bfMeto´do de Euler para h= 0.25’);
legend(’Met-Euler’,’Int Exato’)
grid on
Saı´da
Valor exato da Integral e´: 0.36788
=================================================================
Valor da Integrac¸~ao Numerica de Euler: 0.29630
=================================================================
Erro da Integrac¸~ao Numerica de Euler e´: 19.46 por cento
=================================================================
5
4.2 Rotina me´todo de Euler para h = 0,125
% Exercicio Lista 5 ELIZETE ROCHA
% Me´todo de Euler para integrac¸~ao nume´rica
% Exercicio 2 y’ = -y => y = e^x para h = 0,25
clear all
close all
clc
x0 = 0;
xf = 1;
y0 = 1;
n = 8;
deltaX = (xf - x0)/(n-1);
%y’ = y => y = e^x
x(1) = x0;
y(1) = y0;
for i = 1: n-1
y(i+1) = y(i) + deltaX * -y(i);
x(i+1) = x(i) + deltaX;
end
y_Ex =exp(-x);
er = (y(end) - y_Ex(end)) / y_Ex(end);
Erro =er* (-100);
disp(sprintf(’===========================================================’))
disp(sprintf(’Valor exato da Integral e´: %0.5f ’, y_Ex(end)))
disp(sprintf(’===========================================================’))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Euler: %0.5f ’, y(end)))
disp(sprintf(’===========================================================’))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Euler e´: %0.2f por cento’,Erro))
disp(sprintf(’===========================================================’))
plot(x,y,x,y_Ex,’r-’,’LineWidth’,2)
xlim([0,1]);
title(’\bfMeto´do de Euler para h= 0.125’);
legend(’Met-Euler’,’Int Exato’)
grid on
_________________________________________________________
Saı´da de dados
Valor exato da Integral e´: 0.36788
=================================================================
Valor da Integrac¸~ao Numerica de Euler: 0.33992
=================================================================
Erro da Integrac¸~ao Numerica de Euler e´: 7.60 por cento
6
=================================================================
4.3 Rotina me´todo de Euler para h = 0,0625
% Exercicio Lista 5 ELIZETE ROCHA
% Me´todo de Euler para integrac¸~ao nume´rica
% Exercicio y’ = -y => y = e^x para h = 0,0625
close all
clc
x0 = 0;
xf = 1;
y0 = 1;
n = 16;
deltaX = (xf - x0)/(n-1);
%y’ = y => y = e^x
x(1) = x0;
y(1) = y0;
for i = 1: n-1
y(i+1) = y(i) + deltaX * -y(i);
x(i+1) = x(i) + deltaX;
end
y_Ex =exp(-x);
er = (y(end) - y_Ex(end)) / y_Ex(end);
Erro =er* (-100);
disp(sprintf(’=============================================================’))
disp(sprintf(’Valor exato da Integral e´: %0.5f ’, y_Ex(end)))
disp(sprintf(’============================================================’))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Euler: %0.5f ’, y(end)))
disp(sprintf(’=============================================================’))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Euler e´: %0.2f por cent ’, Erro))
disp(sprintf(’============================================================’))
plot(x,y,x,y_Ex,’r-’,’LineWidth’,2)
xlim([0,1]);
title(’\bfMeto´do de Euler para h= 0.0625’);
legend(’Met-Euler’,’Int. Exata’)
grid on
--------------------------------------------------------------
Saı´da de dados
Valor exato da Integral e´: 0.36788
=================================================================
Valor da Integrac¸~ao Numerica de Euler: 0.35526
=================================================================
7
Erro da Integrac¸~ao Numerica de Euler e´: 3.43 por cento
=================================================================
4.4 Rotina comparac¸a˜o entre o me´todo de Euler e RK4 para
h = 0,25
% Exercicio Lista 5 ELIZETE ROCHA
% Me´todo de Euler e Runge-Kutta 4o ordem para integrac¸~ao nume´rica
% Exercicio 2 y’= -y, y(0)=1 h = 0,25
function metodo_1_EeRK4()
close all
dt = 0.25;
y0 = 1;
% soluc¸~ao analitica
t = 0:dt:1;
y_a_sol = exp(-t);
figure(1)
plot(t,y_a_sol,’k-’,’LineWidth’,2)
hold on
% metodo de EULER
tE(1) = 0;
yE(1) = y0;
n = 1;
while tE < 1
tE(n+1) = tE(n) + dt;
yE(n+1) = yE(n) + f(yE(n),tE(n))*dt;
n = n + 1;
end
figure(1)
plot(tE,yE,’g-’,’LineWidth’,2)
% RK2 opcional no exercico n~ao foi pedido
t2(1) = 0;
y2(1) = y0;
n = 1;
while t2 < 1
t2(n+1) = t2(n) + dt;
k1 = f(y2(n),t2(n));
k2 = f(y2(n)+k1*dt,t2(n) + dt);
phi = 1/2*(k1 + k2);
y2(n+1) = y2(n) + phi*dt;
n = n + 1;
end
8
figure(1)
plot(t2,y2,’r--’,’LineWidth’,2)
% RK4
t4(1) = 0;
y4(1) = y0;
n = 1;
while t4 < 1
t4(n+1) = t4(n) + dt;
k1 = f(y4(n),t4(n));
k2 = f(y4(n)+k1*dt/2,t4(n) + dt/2);
k3 = f(y4(n)+k2*dt/2,t4(n) + dt/2);
k4 = f(y4(n)+k3*dt,t4(n) + dt);
phi = (1/6)*(k1 + 2*k2 + 2*k3 + k4);
y4(n+1) = y4(n) + phi*dt;
n = n + 1;
end
figure(1)
plot(t4,y4,’b-’,’LineWidth’,2)
hold on
legend(’SolucaoAnalitica’,’Euler’,’RK2’,’RK4’)
grid on
title(’\bfMeto´do de Euler para h= 0.25’);
erE = (yE(end) - y_a_sol(end)) / y_a_sol(end);
erRK2 = (y2(end) - y_a_sol(end)) / y_a_sol(end);
erRK4 = (y4(end) - y_a_sol(end)) / y_a_sol(end);
E1=erE *100;
E2 =erRK2*100;
E3 = erRK4*100;
disp(sprintf(’Valor exato da Integral e´: %0.10f ’, y_a_sol(end)))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Euler: %0.10f ’, yE(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Euler e´: %0.10f por cento ’, E1))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Runge-Kutta2: %0.10f ’, y2(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Runge-Kutta2 e´: %0.10f por cento ’, E2))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Runge-Kutta4: %0.10f ’, y4(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Runge-Kutta4 e´: %0.10f por cento ’, E3))
disp(sprintf(’=========================================================================
==============’))
disp(sprintf(’========================= Resumo das respostas para o exerecı´cio 2 =======
=============’))
disp(sprintf(’==========================================================================
=============’))
disp(sprintf(’Valor da Integrac¸~ao Numerica Soluc¸~ao Analitica e´: %0.10f ’, y_a_sol(end)))
disp(sprintf(’===========================================================================
============’))
disp(sprintf(’Erros encontrados; para Euler:%0.4f por cento, e Runger-Kutta4:%0.4f por
cento ’,E1, E3))
9
disp(sprintf(’==========================================================================
=============’))
disp(sprintf(’Valor da Integrac¸~ao Numerica de: Euler e´:%0.10f Runge-Kutta4: %0.10f ’,
yE(end),y4(end)))
disp(sprintf(’=======================================================================
================’))
function ydot = f(y,t)
ydot = -y;
Saı´da dos dados
Valor exato da Integral e´: 0.3678794412
Valor da Integrac¸~ao Numerica de Euler: 0.3164062500
Erro da Integrac¸~ao Numerica de Euler e´: -13.9918640214 por cento
Valor da Integrac¸~ao Numerica de Runge-Kutta2: 0.3725290298
Erro da Integrac¸~ao Numerica de Runge-Kutta2 e´: 1.2638892404 por cento
Valor da Integrac¸~ao Numerica de Runge-Kutta4: 0.3678941994
Erro da Integrac¸~ao Numerica de Runge-Kutta4 e´: 0.0040117043 por cento
=========================================================================
========================= Resumo das respostas para o exerecı´cio 2 =======
==========================================================================
Valor da Integrac¸~ao Numerica Soluc¸~ao Analitica e´: 0.3678794412
=========================================================================
Erros encontrados; para Euler:-13.9919 por cento, e Runger-Kutta4: 0.0040
===========================================================================
Valor da IntNumerica de: Euler e´:0.3164062500 Runge-Kutta4: 0.3678941994
===========================================================================
4.5 Rotina comparac¸a˜o entre o me´todo de Euler e RK4 para
h = 0,0125
% Exercicio Lista 5 ELIZETE ROCHA
% Me´todo de Euler e Runge-Kutta 4oordem para integrac¸~ao nume´rica
% Exercicio 2 y’= -y, y(0)=1, h/2 = 0,125
function metodo_2_EeRK4()
close all
dt = 0.125;
y0 = 1;
% soluc¸~ao analitica
t = 0:dt:1;
y_a_sol = exp(-t);
figure(1)
plot(t,y_a_sol,’k-’,’LineWidth’,2)
hold on
% metodo de EULER
tE(1) = 0;
10
yE(1) = y0;
n = 1;
while tE < 1
tE(n+1) = tE(n) + dt;
yE(n+1) = yE(n) + f(yE(n),tE(n))*dt;
n = n + 1;
end
figure(1)
plot(tE,yE,’g-’,’LineWidth’,2)
% RK2 opcional no exercico n~ao foi pedido
t2(1) = 0;
y2(1) = y0;
n = 1;
while t2 < 1
t2(n+1) = t2(n) + dt;
k1 = f(y2(n),t2(n));
k2 = f(y2(n)+k1*dt,t2(n) + dt);
phi = 1/2*(k1 + k2);
y2(n+1) = y2(n) + phi*dt;
n = n + 1;
end
figure(1)
plot(t2,y2,’r--’,’LineWidth’,2)
% RK4
t4(1) = 0;
y4(1) = y0;
n = 1;
while t4 < 1
t4(n+1) = t4(n) + dt;
k1 = f(y4(n),t4(n));
k2 = f(y4(n)+k1*dt/2,t4(n) + dt/2);
k3 = f(y4(n)+k2*dt/2,t4(n) + dt/2);
k4 = f(y4(n)+k3*dt,t4(n) + dt);
phi = (1/6)*(k1 + 2*k2 + 2*k3 + k4);
y4(n+1) = y4(n) + phi*dt;
n = n + 1;
end
figure(1)
plot(t4,y4,’b-’,’LineWidth’,2)
hold on
legend(’SolucaoAnalitica’,’Euler’,’RK2’,’RK4’)
grid on
title(’\bfMeto´do de Euler para h= 0.125’);
11
erE = (yE(end) - y_a_sol(end)) / y_a_sol(end);
erRK2 = (y2(end) - y_a_sol(end)) / y_a_sol(end);
erRK4 = (y4(end) - y_a_sol(end)) / y_a_sol(end);
E1=erE *100;
E2 =erRK2*100;
E3 = erRK4*100;
disp(sprintf(’Valor exato da Integral e´: %0.10f ’, y_a_sol(end)))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Euler: %0.10f ’, yE(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Euler e´: %0.10f por cento ’, E1))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Runge-Kutta2: %0.10f ’, y2(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Runge-Kutta2 e´: %0.10f por cento
’, E2))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Runge-Kutta4: %0.10f ’, y4(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Runge-Kutta4 e´: %0.10f por cento ’
, E3))
disp(sprintf(’================================================================
=======================’))
disp(sprintf(’========================= Resumo das respostas para o exerecı´cio 2
====================’))
disp(sprintf(’================================================================
======================’))
disp(sprintf(’Valor da Integrac¸~ao Numerica Soluc¸~ao Analitica e´: %0.10f ’,
y_a_sol(end)))
disp(sprintf(’==================================================================
=====================’))
disp(sprintf(’Erros encontrados; para Euler:%0.4f por cento, e Runger-Kutta4: %0.8f
por cento ’,E1, E3))
disp(sprintf(’=====================================================================
==================’))
disp(sprintf(’Valor da Integrac¸~ao Numerica de: Euler e´:%0.10f Runge-Kutta4: %0.10f
’, yE(end),y4(end)))
disp(sprintf(’=====================================================================
==================’))
function ydot = f(y,t)
ydot = -y;
==============================
saı´da de dados
===============================
Valor exato da Integral e´: 0.3678794412
Valor da Integrac¸~ao Numerica de Euler: 0.3436089158
Erro da Integrac¸~ao Numerica de Euler e´: -6.5974128069 por cento
Valor da Integrac¸~ao Numerica de Runge-Kutta2: 0.3689332441
Erro da Integrac¸~ao Numerica de Runge-Kutta2 e´: 0.2864533299 por cento
Valor da Integrac¸~ao Numerica de Runge-Kutta4: 0.3678802719
12
Erro da Integrac¸~ao Numerica de Runge-Kutta4 e´: 0.0002258214 por cento
==============================================================================
========================= Resumo das respostas para o exerecı´cio 2 ===========
==============================================================================
Valor da Integrac¸~ao Numerica Soluc¸~ao Analitica e´: 0.3678794412
==============================================================================
Erros encontrados; para Euler:-6.5974 por cento, e Runger-Kutta4: 0.00022582
==============================================================================
Valor da Integrac¸~ao Numerica de: Euler e´:0.3436089158 Runge-Kutta4: 0.3678802719
==============================================================================
4.6 Rotina comparac¸a˜o entre o me´todo de Euler e RK4 para
h = 0,0625
% Exercicio Lista 5 ELIZETE ROCHA
% Me´todo de Euler e Runge-Kutta 4oordem para integrac¸~ao nume´rica
% Exercicio 2 y’= -y, y(0)=1, h/4 = 0,0625 onede h inicial = 0,25
function metodo_4_EeRK4()
close all
dt = 0.0625;
y0 = 1;
% soluc¸~ao analitica
t = 0:dt:1;
y_a_sol = exp(-t);
figure(1)
plot(t,y_a_sol,’k-’,’LineWidth’,2)
hold on
% metodo de EULER
tE(1) = 0;
yE(1) = y0;
n = 1;
while tE < 1
tE(n+1) = tE(n) + dt;
yE(n+1) = yE(n) + f(yE(n),tE(n))*dt;
n = n + 1;
end
figure(1)
plot(tE,yE,’g-’,’LineWidth’,2)
% RK2 opcional no exercico n~ao foi pedido
t2(1) = 0;
y2(1) = y0;
n = 1;
13
while t2 < 1
t2(n+1) = t2(n) + dt;
k1 = f(y2(n),t2(n));
k2 = f(y2(n)+k1*dt,t2(n) + dt);
phi = 1/2*(k1 + k2);
y2(n+1) = y2(n) + phi*dt;
n = n + 1;
end
figure(1)
plot(t2,y2,’r--’,’LineWidth’,2)
% RK4
t4(1) = 0;
y4(1) = y0;
n = 1;
while t4 < 1
t4(n+1) = t4(n) + dt;
k1 = f(y4(n),t4(n));
k2 = f(y4(n)+k1*dt/2,t4(n) + dt/2);
k3 = f(y4(n)+k2*dt/2,t4(n) + dt/2);
k4 = f(y4(n)+k3*dt,t4(n) + dt);
phi = (1/6)*(k1 + 2*k2 + 2*k3 + k4);
y4(n+1) = y4(n) + phi*dt;
n = n + 1;
end
figure(1)
plot(t4,y4,’b-’,’LineWidth’,2)
hold on
legend(’SolucaoAnalitica’,’Euler’,’RK2’,’RK4’)
grid on
title(’\bfMeto´do de Euler para h= 0.0625’);
erE = (yE(end) - y_a_sol(end)) / y_a_sol(end);
erRK2 = (y2(end) - y_a_sol(end)) / y_a_sol(end);
erRK4 = (y4(end) - y_a_sol(end)) / y_a_sol(end);
E1=erE *100;
E2 =erRK2*100;
E3 = erRK4*100;
disp(sprintf(’Valor exato da Integral e´: %0.10f ’, y_a_sol(end)))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Euler: %0.10f ’, yE(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Euler e´: %0.10f por cento ’, E1))
disp(sprintf(’Valorda Integrac¸~ao Numerica de Runge-Kutta2: %0.10f ’, y2(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Runge-Kutta2 e´: %0.10f por cento ’,
E2))
disp(sprintf(’Valor da Integrac¸~ao Numerica de Runge-Kutta4: %0.10f ’, y4(end)))
disp(sprintf(’Erro da Integrac¸~ao Numerica de Runge-Kutta4 e´: %0.10f por cento ’,
E3))
disp(sprintf(’==================================================================
14
=====================’))
disp(sprintf(’========================= Resumo das respostas para o exerecı´cio 2
====================’))
disp(sprintf(’===================================================================
====================’))
disp(sprintf(’Valor da Integrac¸~ao Numerica Soluc¸~ao Analitica e´: %0.10f ’,
y_a_sol(end)))
disp(sprintf(’=============================================================
==========================’))
disp(sprintf(’Erros encontrados; para Euler:%0.4f por cento, e Runger-Kutta4:
%0.8f por cento ’,E1, E3))
disp(sprintf(’================================================================
=======================’))
disp(sprintf(’Valor da Integrac¸~ao Numerica de: Euler e´:%0.10f Runge-Kutta4: %0.10
f ’, yE(end),y4(end)))
disp(sprintf(’================================================================
=======================’))
function ydot = f(y,t)
ydot = -y;
------------------------------------
Saı´da de dados
===========================
Valor exato da Integral e´: 0.3678794412
Valor da Integrac¸~ao Numerica de Euler: 0.3560741305
Erro da Integrac¸~ao Numerica de Euler e´: -3.2090161609 por cento
Valor da Integrac¸~ao Numerica de Runge-Kutta2: 0.3681305387
Erro da Integrac¸~ao Numerica de Runge-Kutta2 e´: 0.0682553894 por cento
Valor da Integrac¸~ao Numerica de Runge-Kutta4: 0.3678794905
Erro da Integrac¸~ao Numerica de Runge-Kutta4 e´: 0.0000133960 por cento
===============================================================================
========================= Resumo das respostas para o exerecı´cio 2 ===========
===============================================================================
Valor da Integrac¸~ao Numerica Soluc¸~ao Analitica e´: 0.3678794412
================================================================================
Erros encontrados; para Euler:-3.2090 por cento, e Runger-Kutta4: 0.00001340 por
===============================================================================
Valor da Integrac¸~ao Numerica de: Euler e´:0.3560741305 Runge-Kutta4: 0.3678794905
==============================================================================
15