Metodos de Resolução Equações Diferenciais Ordinárias com PVI

Prévia do material em texto

UNIVERSIDADE ESTADUAL DO MARANHÃO
CENTRO DE CIENCIAS TECNOLOGICAS
CURSO DE ENGENHARIA DA COMPUTAÇÃO
Eryck de Araujo Oliveira 11092037
Relatório Técnico Científico 
Métodos Numéricos Para Resolução de EDO 
São Luís – MA
2014
Eryck de Araujo Oliveira 11092037
Relatório Técnico Científico 
Métodos Numéricos Para Resolução de EDO 
Relatório Técnico como requisito parcial para obtenção de aprovação na disciplina Métodos Numéricos Aplicados à Engenharia, no Curso de Engenharia da Computação, na Universidade Estadual do Maranhão.
Prof. Dr. Mariano Costa Amaral
São Luís – MA
2014
INTRODUÇÃO 
Uma equação diferencial com apenas uma variável independente é chamada de equação diferencial ordinária (EDO). Uma EDO de primeira ordem envolve a derivada primeira da variável dependente em relação à variável independente. Para que uma solução específica seja obtida, uma EDO de primeira ordem deve ter uma condição ou restrição inicial que especifique o valor da variável dependente para um valor particular da variável independente. 
Como problemas reais usando EDO para solucioná-los são tipicamente problemas que envolvem a variação temporal (isto é, problemas nos quais o tempo t é a variável independente), a restrição é chamada de condição inicial e o problema é chamado de problema de valor inicial (PVI). 
A solução analítica de uma EDO de primeira ordem é uma expressão matemática da função y(x) que satisfaz a equação diferencial e inclui o valor y(x1) = y1. Uma vez conhecida a função y(x), pode-se calcular o seu valor em qualquer x.
A solução numérica de uma EDO é formada por um conjunto de pontos discretos que representam a função y(x) de forma aproximada. Quando se resolve numericamente uma equação diferencial, o enunciado do problema também inclui o domínio da solução. Daí da solução aproximada de uma EDO clássica que vem os Métodos Numéricos para solução de EDO com PVI e este presente trabalho apresenta os métodos referente às essas soluções bem como as soluções, os erros encontrados em cada um e gráficos mostrando a precisão e a comparação entre os mesmos.
FUNDAMENTAÇÃO TEÓRICA
Método de Euler
 O método explícito de Euler, também chamado de método de Euler progressivo, é uma técnica numérica de passo simples usada na solução de EDOs de primeira ordem. O método de Euler assume que, em uma pequena distância h na vizinhança de (xi, yi), a função y(x) tem uma inclinação constante e igual à inclinação em (xi, yi). A partir dessa hipótese, o próximo ponto da solução numérica (, ) é calculado usando: 
Método de Euler Modificado
O método de Euler modificado é uma técnica numérica explícita de passo simples usada para resolver EDOs de primeira ordem. Ela é uma versão modificada do método explícito de Euler. No método de Euler modificado, a inclinação utilizada no cálculo do valor de é modificada para incluir o efeito de variação desse parâmetro ao longo do intervalo. A inclinação usada no método de Euler modificado é a média da inclinação no início do intervalo e de uma estimativa para a inclinação no final do intervalo.
Uma vez calculadas as duas inclinações, uma melhor estimativa para o valor de é calculada a partir da média dessas inclinações:
Método de Runge-Kutta
O método de Runge-Kutta é provavelmente um dos métodos mais populares. O método de Runge-Kutta de quarta ordem também é um dos mais preciosos para obter soluções aproximadas de valor inicial. Cada método de Runge-Kutta consiste em comparar um polinômio de Taylor apropriado para eliminar o cálculo das derivadas. Fazendo-se várias avaliações da função a cada passo. Estes métodos podem ser construídos para qualquer ordem α. O qual será abordado nesse trabalho. O método tem a forma da expressão abaixo:
Métodos Multipassos (Adams-Bashforth)
Em métodos de passo simples, a solução no ponto é obtida usando apenas o valor de e no ponto anterior. Em métodos multipasso, a solução no ponto é calculada a partir de dois ou mais pontos anteriores. Métodos multipasso podem ser implícitos ou explícitos. Em métodos multipasso explícitos, a solução é calculada a partir de uma fórmula explícita. Normalmente na definição de dois ou mais pontos anteriores são utilizados métodos de passo simples para a determinação dos mesmos. As fórmulas do método de Adams-Bashfort são deduzidas a partir da integração da equação diferencial ao longo de um intervalo arbitrário [xi, xi + 1]: 
Método das Diferenças Finitas
 Este método consiste em escrever a EDO dada e suas condições de contorno na forma de equações em diferenças finitas e resolver o sistema de equações algébricas resultante. O domínio da solução [a, b] é dividido em N subintervalos de mesma largura h, que são definidos por (N + 1) pontos chamados de pontos de malha (em geral, subintervalos podem ter larguras desiguais).
A largura de cada subintervalo é portanto h = (b – a)/N. Os pontos a e b são os pontos finais, e o restante dos pontos são os pontos internos. A equação diferencial é então escrita em cada um dos pontos internos do domínio.
METODOLOGIA
Neste presente trabalho foi usada a ferramenta computacional Matlab, formulando os códigos dos métodos a partir da fundamentação teórica obtida dos mesmos. Após essa etapa começou-se a etapa de testes e comparação dos métodos, onde foram escolhidas três Equações Diferenciais Ordinárias arbitrárias com condições iniciais, adotando-se um valor inicial e final para o parâmetro x e um valor para a condição inicial y. 
Então os códigos foram formulados todos com a mesma métrica obedecendo o padrão function [x,y] = nome_do_método(f,xini,xfinal,yini,n) onde f = função y'=f(x,y), yini = condição inicial, xini = limite inferior, xfinal = limite superior, n = número de passos desejado.
E com a função denominada erro se obtém a precisão de cada método em relação à solução analítica da EDO utilizada nos testes. Assim pode-se obter a comparação entre os métodos e ver qual é o mais preciso, qual obtém o menor custo computacional e o menor tempo. 
RESULTADOS 
Com a função com x(0)=0, x=2.5, y(0)=3 e h=0.1. Com solução exata em .
Temos os seguintes resultados processados pela ferramenta Matlab: 
	Iterações
	Analítica
	Euler
	Euler_mod
	Runge-kutta
	Multipassos
	
	3.0000
	3.0000
	3.0000
	3.0000
	3.0000
	
	3.3104
	3.3400
	3.3093
	3.3104
	3.3104
	
	3.5665
	3.6185
	3.5645
	3.5665
	3.5665
	
	3.7750
	3.8435
	3.7724
	3.7750
	3.7750
	
	3.9419
	4.0221
	3.9387
	3.9419
	3.9419
	
	4.0723
	4.1603
	4.0689
	4.0723
	4.0723
	
	4.1710
	4.2635
	4.1673
	4.1709
	4.1710
	
	4.2419
	4.3366
	4.2382
	4.2419
	4.2419
	
	4.2888
	4.3836
	4.2851
	4.2888
	4.2888
	
	4.3149
	4.4082
	4.3112
	4.3149
	4.3149
	
	4.3229
	4.4136
	4.3193
	4.3229
	4.3229
	
	4.3153
	4.4025
	4.3119
	4.3153
	4.3153
	
	4.2944
	4.3775
	4.2911
	4.2944
	4.2944
	
	4.2620
	4.3406
	4.2590
	4.2620
	4.2620
	
	4.2199
	4.2936
	4.2170
	4.2199
	4.2199
	
	4.1696
	4.2383
	4.1669
	4.1696
	4.1696
	
	4.1123
	4.1761
	4.1098
	4.1123
	4.1123
	
	4.0493
	4.1081
	4.0470
	4.0493
	4.0493
	
	3.9815
	4.0355
	3.9794
	3.0000
	3.9815
	
	3.4361
	3.4619
	3.4352
	3.4361
	3.4361
Figura 1: ilustração da curva da solução analítica (função 1) com as soluções numéricas, com seus erros relativos.
Para a função com x(0)=0, x=1, y(0)=2 e h=0.01. Com solução exata em .
Temos os seguintes resultados processados pela ferramenta Matlab: 
	Iterações
	Analítica
	Euler
	Euler_mod
	Runge-kutta
	Multipassos
	1.
	2.0000
	2.0000
	2.0000
	2.0000
	2.0000
	10.
	2.2950
	2.2921
	2.2950
	2.2950
	2.2950
	20.
	2.6853
	2.6785
	2.6853
	2.6853
	2.6853
	30.
	3.1476
	3.1362
	3.1476
	3.1476
	3.1476
	40.
	3.6888
	3.6718
	3.6887
	3.6888
	3.6888
	50.
	4.3160
	4.2925
	4.3159
	4.3160
	4.3160
	60.
	5.0370
	5.00605.0369
	5.0370
	5.0370
	70.
	5.8600
	5.8203
	4.2382
	5.8600
	5.8600
	80.
	6.7940
	6.7440
	6.7938
	6.7940
	6.7940
	90.
	7.8484
	7.7866
	7.8481
	7.8484
	7.8484
	100.
	9.1596
	9.0828
	9.1593
	9.1596
	9.1596
Figura 2: ilustração da curva da solução analítica (função 2) com as soluções numéricas, com seus erros relativos.
Para a função com x(0)=0, x=1, y(0)=0 e h=0.05. Com solução exata em .
Temos os seguintes resultados processados pela ferramenta Matlab: 
	Iterações
	Analítica
	Euler
	Euler_mod
	Runge-kutta
	Multipassos
	
	0
	0
	0
	0
	0
	
	0.0053
	0
	 0.0063
	0.0053
	0.0053
	
	0.0184
	0.0125
	 0.0195
	0.0184
	0.0184
	
	0.0362
	0.0313
	 0.0372
	0.0362
	0.0362
	
	0.0568
	0.0531
	 0.0576
	0.0568
	0.0567
	
	0.0791
	0.0766
	 0.0798
	0.0791
	0.0790
	
	0.1025
	0.1008
	 0.1030
	0.1025
	0.1024
	
	0.1265
	0.1254
	 0.1269
	0.1265
	0.1265
	
	0.1509
	0.1502
	 0.1512
	0.1509
	0.1509
	
	0.1756
	0.1751
	 0.1757
	0.1756
	0.1755
	
	0.2003
	0.2000
	 0.2005
	0.2003
	0.2003
	
	0.2252
	0.2250
	 0.2253
	0.2252
	0.2252
	
	0.2501
	0.2500
	 0.2502
	0.2501
	0.2501
	
	0.2751
	0.2750
	 0.2751
	0.2751
	0.2751
	
	0.3000
	0.3000
	 0.3001
	0.3000
	0.3000
	
	0.3250
	0.3250
	 0.3250
	0.3250
	0.3250
	
	0.3500
	0.3500
	 0.3500
	0.3500
	0.3500
	
	0.3750
	0.3750
	 0.3750
	0.3750
	0.3750
	
	0.4000
	0.4000
	 0.4000
	0.4000
	0.4000
	
	0.4250
	0.4250
	 0.4250
	0.4250
	0.4250
Figura 2: ilustração da curva da solução analítica (função 3) com as soluções numéricas, com seus erros relativos.
CONCLUSÃO
Com base nos resultados obtidos vemos que cada função nos dá uma solução analítica e uma solução numérica obtida por meio de cada um dos métodos. Depois de processadas as soluções viu-se que o método de Euler nos dá uma solução aproximada mas distante da solução exata, já o método de Euler Modificado se aproxima mais da solução exata com o mesmo número de iterações. Com o mesmo número de iterações os métodos de Multipassos e de Runge-Kutta nos deram a maior precisão com erros na casa de na primeira função, na segunda função e na terceira função. 
Nos dando base da sua superioridade em relação aos métodos explícitos de Euler e Euler Modificado. Sendo eles a base das soluções de EDOs dos mais variados tipos.
ANEXOS
Código Fonte (Euler Explicito) 
function [x,y]=euler(f,xini,yini,xfinal,h)
% Aproximação de Euler para o problema de valor inicial de uma EDO
% Método Explicito de Euler
% Calculo do h de xini, xfinal e n
n=(xfinal-xini)/h;
% Initialização de x e y como vector de colunas
x=[xini zeros(1,n)];
y=[yini zeros(1,n)];
% Calculo de x e y
for i=1:n
 x(i+1)=x(i)+h;
 y(i+1)=y(i)+h*f(x(i),y(i));
 end
end
Código Fonte (Euler Implícito) 
function [x,y]=euler_imp(f,xini,yini,xfinal,h)
n=(xfinal-xini)/h;
% Initialização de x e y como vector de colunas
x=[xini zeros(1,n)];
y=[yini zeros(1,n)];
% Calculo de x e y
for i=1:n
 x(i+1)=x(i)+h;
 ynew=y(i)+h*(f(x(i),y(i)));
 y(i+1)=y(i)+h*f(x(i+1),ynew);
end
end
Código Fonte (Euler Modificado) 
function [x,y]=euler_mod(f,xini,yini,xfinal,h)
% Aproximação de Euler para o problema de valor inicial de uma EDO
% Método de Euler Modificado
% Calculo do h de xini, xfinal e n
n=(xfinal-xini)/h;
% Initialização de x e y como vector de colunas
x=[xini zeros(1,n)];
y=[yini zeros(1,n)];
% Calculo de x e y
for i=1:n
 x(i+1)=x(i)+h;
 ynew=y(i)+h*f(x(i),y(i));
 y(i+1)=y(i)+(h/2)*(f(x(i),y(i))+f(x(i+1),ynew));
 end
end
Código Fonte (Runge-Kutta 4ª Ordem) 
function [x,y] = RK4(f,xini,yini,xfinal,h)
% edoRK4 resolve uma EDO de primeira ordem com valor inicial usando o
% método de Runge-Kutta de quarta ordem.
n = (xfinal - xini) / h;
halfh = h / 2;
y(1,:) = yini;
x(1) = xini;
h6 = h/6;
for i = 1 : n
 x(i+1) = x(i) + h;
 th2 = x(i) + halfh;
 K1 = f(x(i), y(i,:));
 K2 = f(th2, y(i,:) + halfh * K1);
 K3 = f(th2, y(i,:) + halfh * K2);
 K4 = f(x(i+1), y(i,:) + h * K3);
 y(i+1,:) = y(i,:) + (K1 + K2+K2 + K3+K3 + K4) * h6;
end;
Código Fonte (Multipasso Adams Bashforth Moulton) 
function [x,y] = Multi_abm(f,xini,yini,xfinal,h)
%Processo preditor corrector de Adams-Bashforth-Moulton (de quarta ordem)
n=((xfinal-xini)/h);
% Definir vectores coluna de y e x
x=zeros(1,n+1);
y=zeros(1,n+1);
y(1)=yini;
x=xini:h:xfinal;
%cálculo de y
%As 4 primeiras coordenadas de x e y devem ser obtidas por RK4
for i=1:3
 F1=feval(f,x(i),y(i));
 F2=feval(f,x(i)+(h/2),y(i)+(h*F1/2));
 F3=feval(f,x(i)+(h/2),y(i)+(h*F2/2));
 F4=feval(f,x(i)+h, y(i)+h*F3);
 y(i+1)=y(i)+(h/6)*(F1+2*F2+2*F3+F4);
 for j=4:n
 %função do preditor
 G1=feval(f,x(j),y(j));
 G2=feval(f,x(j-1),y(j-1));
 G3=feval(f,x(j-2),y(j-2));
 G4=feval(f,x(j-3),y(j-3));
 p=y(j)+(h/24)*(55*G1-59*G2+37*G3-9*G4);
 %função do corretor
 y(j+1)=y(j)+(h/24)*(9*feval(f,x(j+1),p)+19*G1-5*G2+G3);
 end
end
 REFERÊNCIAS
[1] PRESS, W. H.; TEUKOLSKY, S. A.; VETTERLING W. T.; FLANNERY, B. P., “Numerical Recipes in C: The Art of Scientific Computing”. Cambridge University Press, 1992. 
[2] GILAT, Amos Pozo. ”Métodos numéricos para engenheiros e cientistas [recurso eletrônico] : uma introdução com aplicações usando o MATLAB”; tradução Alberto Resende de Conti. – Dados eletrônicos. – Porto Alegre: Bookman, 2008.
[3] DÉCIO, Esperandio; MENDES, João T.; MONKEN, Luiz Henry,” Cálculo Numérico”, Pearson/Pretice Hall, 2005.
[4] BOYCE W. E.; PRIMA, R. C. Di. “Elementary Di®erential Equations and Boundary Value Problems”. John Wiley and Sons, 2001.
[5] COOMBES, K. R., HUNT, B. R., LIPSMAN, R. L. “Diferential Equations with MATLAB”. John Wiley and Sons, 2000.
[6] NAKAMURA, S. “Numerical Analysis with MATLAB”. Prentice Hall, 2002.
[7] MOLER, C. B. “Numerical Computing with MATLAB”. Siam, 2004.