Trabalho Runge Kutta (com programas)

Prévia do material em texto

Trabalho 1
Luan Bottin De Toni
IF-UFRGS
30 de junho de 2017
1 Introdução
Este trabalho tem como objetivo analisar brevemente os métodos de
Runge-Kutta para a resolução numérica de soluções de equações diferenciais
ordinárias, estudados em aula. Será comparado os resultados do método de
Runge-Kuta de segunda e quarta ordens, juntamente com o método de Euler-
Cromer estudado anteriormente, para um problema com forças dependentes
da velocidade.
2 Lançameto Oblíquo
Considere uma bola de massam que foi arremessada com uma velocidade
inicial ~v0 = 60m/s e ângulo θ0 = 60
◦
. Esta bola está sujeita à ação da força
gravitacional e força de arraste. Ou seja:
~Fres = −m~g −mγ~v (1)
onde γ é a constante de amortecimento.
O sistema a ser resolvido para aplicação do algoritmo é:
x˙ = vx (2)
y˙ = vy (3)
v˙x = −γx˙ (4)
v˙y = −g − γy˙ (5)
A fim de verificar a qualidade dos resultados numéricos obtidos, será
comparado o resultado utilizando os métodos numéricos com o resultado
analítico para a posição x(t) e y(t), expressados por:
x(t) =
V0x
γ
(1− e−γt) (6)
y(t) = −gt
γ
− (γV0y + g
γ2
)(1− e−γt) (7)
1
A fim de analisar o efeito do parâmetro h nos algoritmos, buscou-se variar
tal parâmetro para os métodos RK2, RK4 e Euler-Cromer para γ = 0.04.
Calculando o erro relativo dos métodos para x(t) com auxílio da equação (4)
plotou-se um gráfico de h x Erro relativo (código do programa na seção 4.1),
aplicando log nos dois eixos para melhor visualização:
 1e−08
 1e−07
 1e−06
 1e−05
 0.0001
 0.001
 0.01
 0.1
 1e−05 0.0001 0.001 0.01 0.1 1
Er
ro
 re
la
tiv
o
h
Erro x h 
γ=0.04 v0=600m/s θ0=60°
RK4
RK2
Euler Cromer
Figura 1: Erro relativo para diferentes valores de h.
Os valores gerados correspondem ao momento que a bola atinge o chão.
Ao gerar os valores para plotar o gráfico acima, notou-se que o tempo de
processamento é relativamente demorado (na ordem de segundos), isto ocorre
pois os passos t feitos pelo programa variam da forma t = t + h, quando h
é muito pequeno é necessário uma enorme quantidade de processos. Nota-
se, também, que para γ = 0.04 os métodos RK2 e RK4 apresentam erros
menores.
Para analisar o efeito do parâmetro γ com o método RK4, será escolhido
o parâmetro h = 0.01 tendo em vista o número pequeno de processos que
este exige e seu erro relativo aceitável (da ordem de 10−5).
O gráfico da posição em x pela posição em y (código do programa utili-
zado na seção 4.2), para diferentes γ encontra-se a seguir:
2
 0
 2000
 4000
 6000
 8000
 10000
 12000
 0 5000 10000 15000 20000 25000
y(t
)
x(t)
X x Y (RUNGE-KUTTA 4)
γ=0.005
γ=0.01
γ=0.02
γ=0.04
γ=0.08
Figura 2: Gráfico da posição horizontal e vertical da bola.
Nota-se que quanto maior o γ menor é o alcance da bola, tanto em x
quanto em y, pois a força de arrasto será maior.
2.1 Problema do Scherer
O livro "Física Computacional", do Scherer, trata de um problema seme-
lhante, porém as forças neste caso são dependentes do módulo da velocidade
v2 = v2x + v
2
y . Logo:
x˙ = vx (8)
y˙ = vy (9)
v˙x = −γvx˙ (10)
v˙y = −mg − γvy˙ (11)
Com os mesmos parâmetros utilizados anteriormente, em = 1kg, plotou-
se um gráfico de x(t) por y(t) para o problema do Scherer.
3
 0
 50
 100
 150
 200
 250
 300
 350
 400
 450
 0 50 100 150 200 250 300 350 400 450
y(t
)
x(t)
X x Y (RUNGE-KUTTA 4 -> Scherer)
γ=0.005
γ=0.01
γ=0.02
γ=0.04
γ=0.08
Figura 3: Gráfico da posição horizontal e vertical da bola para o problema
do Scherer.
Como as forças dependem também do módulo da velocidade os alcances
são bem menores que o exemplo anterior. Também vemos uma "queda"mais
aguda que no exemplo anterior (Figura 2) onde forma quase uma parábola.
Pensando em uma situação real, o gráfico do problema do Scherer faz mais
sentido, pois no começo de sua trajetória, a velocidade horizontal está mais
próxima da velocidade inicial; e quando a bola começa sua queda sua velo-
cidade em x já está menor devido às forças de arrasto, isto explica o alcance
maior na primeira metade da trajetória e a queda abrupta da bola no final.
3 Questões
3.1 Como avaliar o erro cometido quando não sabemos a
solução exata?
Para plotar o gráfico da figura 3 utilizou-se o mesmo h da figura 2, por
não termos em mãos a solução analítica do problema tornando impossível a
análise de erro feita anteriormente utilizando o valor analítico.
É possível estimar o erro relativo do passo em execução mesmo sem a
solução analítica. Uma maneira de fazer isso é comparar o valor calculado
do método usado com o valor de um método superior, sendo essa diferença
o erro corrente. Outra maneira é calcular o valor x(t+ h) a partir de x(t) e,
além disso, calcular novamente em dois passos, cada um com h/2; a diferença
4
desses dois valores estima o erro de truncamento local.
Dessa maneira podemos impor um erro admissível no algoritmo e adaptar
o método para que o programa decida qual h usar em cada passo a partir do
erro calculado com uma das maneiras descritas acima.
3.2 O que deveria ser modificado para que as colisões com o
chão fossem levadas em consideração?
O programa feito aqui calcula a trajetória da bola até que a altura chegue
próximo de zero. Porém, seria possível que ele continuasse a calcular as
posições considerando uma colisão elástica com o chão.
Neste caso, no momento em que a bola colide com o chão (y ' 0) sua
velocidade vertical troca de sinal e o programa calcula novamente as posições
plotando várias trajetórias, tendo como condições iniciais as condições finais
da trajetória anterior, com a velocidade em y multiplicada por −1. Porém,
é necessário que o algoritmo possua uma condição de parada como, por
exemplo, um alcance mínimo em x, definido arbitrariarmente; caso contrário
este pode ficar calculando as posições infinitamente. Este alcance mínimo
pode ser definido como, por exemplo, um décimo do alcance inicial (primeira
trajetória da bola) para ter uma visão razoável do gráfico.
4 Códigos
4.1 Programa para análise do erro
#include <stdio.h>
#include <math.h>
double _gamma=0.04, _g=10.0, _m=1.0;
void fRK2(double theta0, double x0, double y0, double v0);
void fRK4(double theta0, double x0, double y0, double v0);
void fEulerCromer(double theta0, double x0, double y0, double v0);
double fAcelx(double vx, double vy);
double fAcely(double vx, double vy);
double fAnaliticox(double t, double theta0, double v0);
double fAnaliticoy(double t, double theta0, double v0);
main () {
double theta0=M_PI/3.0, v0=600.0;
double x0=0.0, y0=0.0;
fRK2(theta0,x0,y0,v0);
5
fRK4(theta0,x0,y0,v0);
fEulerCromer(theta0,x0,y0,v0);
}
//RUNGE-KUTTA 2
void fRK2(double theta0, double x0, double y0, double v0)
{
double t=0,t2=0,h,ax2,ay2;
double vx, vy, vx2, vy2;
double x,y,x2,y2;
FILE *p;
p=fopen("RK2_Erro.dat","w");
for(h=pow(10,-5);h<=1.01;h*=2)
{
vx = v0*cos(theta0);
vy = v0*sin(theta0);
x=x0;
y=y0;
t=0;
while (y>-0.01)
{
t2 = t + h*0.5;
x2 = x + vx*h/2.0;
vx2= vx + fAcelx(vx, vy)*h/2.0;
y2 = y + vy*h/2.0;
vy2= vy + fAcely(vx, vy)*h/2.0;
t+=h;
x = x + vx2*h;
vx= vx + fAcelx(vx2, vy2)*h;
y = y + vy2*h;
vy= vy + fAcely(vx2, vy2)*h;
}
fprintf(p,"%.12e %.12e %.12e %.12e %.12e\n",t,h,x,y,fabs(x-fAnaliticox(t-h,theta0,v0))/fabs(fAnaliticox(t-h,theta0,v0)));
}
6
fclose(p);
}
//RUNGE-KUTTA 4
void fRK4(double theta0, double x0, double y0, double v0)
{
double t=0,h;
double x,y,vx,vy;
double xm,vmx,vmy;
double kx[5], kvx[5], ky[5], kvy[5];
FILE *p;
p=fopen("RK4_Erro.dat","w");
for(h=pow(10,-5);h<=1.01;h*=2)
{
x=x0;
y=y0;
vx=v0*cos(theta0);
vy=v0*sin(theta0);t=0;
while (y>-0.01)
{
//Bloco 1
// X
kx[1]=h*vx;
kvx[1]=h*fAcelx(vx, vy);
// Y
ky[1]=h*vy;
kvy[1]=h*fAcely(vx, vy);
//Bloco 2
// X
xm = x + 0.5*kx[1];
vmx = vx + 0.5*kvx[1];
kx[2] = h*vmx;
// Y
xm = y + 0.5*ky[1];
7
vmy = vy + 0.5*kvy[1];
ky[2] = h*vmy;
// velocidades
kvx[2] = h*fAcelx(vmx, vmy);
kvy[2] = h*fAcely(vmx, vmy);
//Bloco 3
// X
xm = x + 0.5*kx[2];
vmx = vx + 0.5*kvx[2];
kx[3] = h*vmx;
// Y
xm = y + 0.5*ky[2];
vmy = vy + 0.5*kvy[2];
ky[3] = h*vmy;
// velocidades
kvx[3] = h*fAcelx(vmx,vmy);
kvy[3] = h*fAcely(vmx,vmy);
//Bloco 4
// X
xm = x + kx[3];
vmx = vx + kvx[3];
kx[4] = h*vmx;
// Y
xm = y + ky[3];
vmy = vy + kvy[3];
ky[4] = h*vmy;
// velocidades
kvx[4] = h*fAcelx(vmx, vmy);
kvy[4] = h*fAcely(vmx,vmy);
//Finalmente
// X
x = x + (kx[1] + 2*kx[2] + 2*kx[3] + kx[4])/6.0;
vx = vx + (kvx[1] + 2*kvx[2] + 2*kvx[3] + kvx[4])/6.0;
// Y
y = y + (ky[1] + 2*ky[2] + 2*ky[3] + ky[4])/6.0;
vy = vy + (kvy[1] + 2*kvy[2] + 2*kvy[3] + kvy[4])/6.0;
t = t + h;
}
fprintf(p,"%.12e %.12e %.12e %.12e %.12e\n",t,h,x,y,fabs(x-fAnaliticox(t-h,theta0,v0))/fabs(fAnaliticox(t-h,theta0,v0)));
}
8
fclose(p);
}
void fEulerCromer(double theta0, double x0, double y0, double v0)
{
double x,y,vx,vy,t,h;
FILE *p;
p=fopen("EulerCromer_Erro.dat","w");
for(h=pow(10,-5);h<=1.01;h*=2)
{
x=x0;
y=y0;
vx=v0*cos(theta0);
vy=v0*sin(theta0);
t=0;
while(y>-0.01)
{
vx = vx + fAcelx(vx,vy)*h;
x = x + vx*h;
vy = vy + fAcely(vx,vy)*h;
y = y + vy*h;
t+=h;
}
fprintf(p,"%.12e %.12e %.12e %.12e %.12e\n",t,h,x,y,fabs(x-fAnaliticox(t-h,theta0,v0))/fabs(fAnaliticox(t-h,theta0,v0)));
}
fclose(p);
}
double fAcelx(double vx, double vy)
{
double v=1.0;
//v = sqrt(vx*vx + vy*vy); //scherer
return (-_gamma*v*vx);
}
double fAcely(double vx, double vy)
9
{
double v=1.0;
//v = sqrt(vx*vx + vy*vy); //scherer
return(-_g*_m - _gamma*v*vy);
}
double fAnaliticox(double t, double theta0, double v0)
{
return (v0*cos(theta0)/_gamma)*(1-exp(-t*_gamma));
}
double fAnaliticoy(double t, double theta0, double v0)
{
return ( -_g*t/_gamma + ( (_gamma*v0*sin(theta0) + _g)/(_gamma*_gamma) )*(1-exp(-t*_gamma)) );
}
4.2 Programa do lançamento oblíquo com método RK4
#include <stdio.h>
#include <math.h>
double _h=0.01, _gamma=0.08, _g=10.0, _m=1.0;
void fRK4(double theta0, double x0, double y0, double v0);
double fAcelx(double vx, double vy);
double fAcely(double vx, double vy);
main () {
double theta0=M_PI/3.0, v0=600.0;
double x0=0.0, y0=0.0;
fRK4(theta0,x0,y0,v0);
}
//RUNGE-KUTTA 4
void fRK4(double theta0, double x0, double y0, double v0)
{
double t=0;
double x,y,vx,vy;
10
double xm,vmx,vmy;
double kx[5], kvx[5], ky[5], kvy[5];
FILE *p;
p=fopen("RK4Scherer_0.08.dat","w");
x=x0;
y=y0;
vx=v0*cos(theta0);
vy=v0*sin(theta0);
t=0;
while (y>-0.01) //até enconstar no chão
{
fprintf(p,"%.12e %.12e %.12e\n",t,x,y);
//Bloco 1
// X
kx[1]=_h*vx;
kvx[1]=_h*fAcelx(vx, vy);
// Y
ky[1]=_h*vy;
kvy[1]=_h*fAcely(vx, vy);
//Bloco 2
// X
xm = x + 0.5*kx[1];
vmx = vx + 0.5*kvx[1];
kx[2] = _h*vmx;
// Y
xm = y + 0.5*ky[1];
vmy = vy + 0.5*kvy[1];
ky[2] = _h*vmy;
// velocidades
kvx[2] = _h*fAcelx(vmx, vmy);
kvy[2] = _h*fAcely(vmx, vmy);
//Bloco 3
// X
xm = x + 0.5*kx[2];
vmx = vx + 0.5*kvx[2];
kx[3] = _h*vmx;
// Y
11
xm = y + 0.5*ky[2];
vmy = vy + 0.5*kvy[2];
ky[3] = _h*vmy;
// velocidades
kvx[3] = _h*fAcelx(vmx,vmy);
kvy[3] = _h*fAcely(vmx,vmy);
//Bloco 4
// X
xm = x + kx[3];
vmx = vx + kvx[3];
kx[4] = _h*vmx;
// Y
xm = y + ky[3];
vmy = vy + kvy[3];
ky[4] = _h*vmy;
// velocidades
kvx[4] = _h*fAcelx(vmx, vmy);
kvy[4] = _h*fAcely(vmx,vmy);
//Finalmente
// X
x = x + (kx[1] + 2*kx[2] + 2*kx[3] + kx[4])/6.0;
vx = vx + (kvx[1] + 2*kvx[2] + 2*kvx[3] + kvx[4])/6.0;
// Y
y = y + (ky[1] + 2*ky[2] + 2*ky[3] + ky[4])/6.0;
vy = vy + (kvy[1] + 2*kvy[2] + 2*kvy[3] + kvy[4])/6.0;
t = t + _h;
}
fclose(p);
}
double fAcelx(double vx, double vy)
{
double v=1.0;
v = sqrt(vx*vx + vy*vy); //scherer
return (-_gamma*v*vx);
}
12
double fAcely(double vx, double vy)
{
double v=1.0;
v = sqrt(vx*vx + vy*vy); //scherer
return(-_g*_m - _gamma*v*vy);
}
13