Projétil com amortecimento (RK4)

Prévia do material em texto

Me´todos Computacionais da F´ısica
Trabalho 1
Fabio Rasera Figueiredo
244430
25 de marc¸o de 2016
1 Exerc´ıcio 1
Aqui se resolve uma equac¸a˜o diferencial que descreve a trajeto´ria de um
proje´til lanc¸ado obliquamente sob efeito das forc¸as gravitacional e de arrasto
do ar. As forc¸as agindo sa˜o:
~Fres = −m~g −mγ~v (1)
Onde a resoluc¸a˜o se simplifica a resolver as duas equac¸o˜es diferenciais:
x¨ = γx˙ (2)
y¨ = −g − γy˙ (3)
Abaixo o co´digo contendo os me´todos nume´ricos de Euler-Cromer, RK2
(Ponto me´dio) e RK4 (me´todo de Runge) aplicados na resoluc¸a˜o do pro-
blema:
//NOME DO ARQUIVO: projetil.c
//AUTOR: Fabio Rasera Figueiredo
#include <stdio.h>
#include <math.h>
double acc(double gamma1, double grav, double pos, double vel);
double fx(double gamma2, double t1);
double fy(double gamma3, double t2);
int main(void){
double t, E, E1, E0, theta, a0, a, h, x0, r[2], x, y, v0, y0, v[2], vx0, vy0, vx, vy, g[2];
double m1[2], m2[2], m3[2], m4[2], k1[2], k2[2], k3[2], k4[2], gamma;
//Dados iniciais
theta=(M_PI/3.0);
v0=600.0;
h=0.01;
gamma=0.1;
x0=0.0;
y0=0.0;
//Componentes em x e y da velocidade
vx0=v0*cos(theta);
vy0=v0*sin(theta);
1
//-----------RK2----------//
FILE *p;
p = fopen("RK2arrasto.txt","w");
//Dados iniciais
g[0]=0.0; g[1]=10.0; // vetor acelerac¸~ao gravitacional
r[0]=x0; r[1]=y0; // vetor posic¸~ao
v[0]=vx0; v[1]=vy0; // vetor velocidade
t=0; //tempo incial
do{
fprintf(p,"%.12e %.12e %.12e %.12e %.12e %.12e %.12e %.12e\n", r[0],r[1],v[0],v[1],h,t,fx(gamma,t),fy(gamma,t));
//Passo 1 // k1=f(t,x);
{
k1[0]= v[0];
k1[1]= v[1];
m1[0]= acc(gamma, g[0], r[0], k1[0]);
m1[1]= acc(gamma, g[1], r[1], k1[1]);
}
//Passo 2 // k2=f(t+h , x+k1*h*0.5);
{
k2[0]= v[0] + m1[0]*h*0.5;
k2[1]= v[1] + m1[1]*h*0.5;
m2[0]= acc(gamma, g[0], r[0] + k1[0]*h*0.5, k2[0]);
m2[1]= acc(gamma, g[1], r[1] + k1[1]*h*0.5, k2[1]);
}
r[0]= r[0] + k2[0]*h;
r[1]= r[1] + k2[1]*h;
v[0]= v[0] + m2[0]*h;
v[1]= v[1] + m2[1]*h;
t+=h;
}while(r[1]>=0); //Repetir enquanto altura do proje´til for maior que zero
fclose(p);
//-----------RK4----------//
FILE *p1;
p1 = fopen("RK4arrasto.txt","w");
//Reinicializac¸~ao dos dados
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
t=0;
do{
fprintf(p1,"%.12e %.12e %.12e %.12e %.12e %.12e %.12e %.12e\n", r[0],r[1],v[0],v[1],h,t,fx(gamma,t),fy(gamma,t));
//Passo 1 // k1=f(t,x);
{
k1[0]= v[0];
k1[1]= v[1];
m1[0]= acc(gamma, g[0], r[0], k1[0]);
m1[1]= acc(gamma, g[1], r[1], k1[1]);
}
//Passo 2 // k2=f(t+0.5*h , x+0.5*k1*h);
{
k2[0]= v[0] + m1[0]*0.5*h;
k2[1]= v[1] + m1[1]*0.5*h;
m2[0]= acc(gamma, g[0], r[0]+0.5*k1[0]*h, k2[0]);
m2[1]= acc(gamma, g[1], r[1]+0.5*k1[1]*h, k2[1]);
}
//Passo 3 // k3=f(t+0.5*h, x+0.5*k2*h);
{
k3[0]= v[0] +m2[0]*0.5*h;
k3[1]= v[1] +m2[1]*0.5*h;
m3[0]= acc(gamma, g[0], r[0]+0.5*k2[0]*h, k3[0]);
m3[1]= acc(gamma, g[1], r[1]+0.5*k2[1]*h, k3[1]);
}
//Passo 4 // k4=f(t+h, x+k3*h);
{
k4[0]= v[0] +m3[0]*h;
k4[1]= v[1] +m3[1]*h;
2
m4[0]= acc(gamma, g[0], r[0]+k3[0]*h, k4[0]);
m4[1]= acc(gamma, g[1], r[1]+k3[1]*h, k4[1]);
}
//Completando o passo
{
r[0] = r[0] + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])*h/6.0;
r[1] = r[1] + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])*h/6.0;
v[0] = v[0] + (m1[0] + 2*m2[0] + 2*m3[0] + m4[0])*h/6.0;
v[1] = v[1] + (m1[1] + 2*m2[1] + 2*m3[1] + m4[1])*h/6.0;
t+=h;
}
}while(r[1]>=0); //Repetir enquanto a altura dp proje´til for maior ou igual a` zero
fclose(p1);
//------Euler-Cromer------//
FILE *p2;
p2=fopen("EulerCromer.txt","w");
//Reinicializac¸~ao dos dados
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
t=0;
do{
fprintf(p2,"%.12e %.12e %.12e %.12e %.12e %.12e %.12e %.12e\n", r[0],r[1],v[0],v[1],h,t,fx(gamma,t),fy(gamma,t));
v[0]+=acc(gamma, g[0], v[0], v[0])*h;
v[1]+=acc(gamma, g[1], v[1], v[1])*h;
r[0]+=v[0]*h;
r[1]+=v[1]*h;
t+=h;
}while(r[1]>=0); //Repetir enquanto a altura dp proje´til for maior ou igual a` zero
fclose(p2);
return 0;
}
//func¸~ao acelerac¸~ao
double acc(double gamma1, double grav, double pos, double vel){
double a = -grav -gamma1*vel;
return(a);
}
//func¸~ao exata em x
double fx(double gamma2, double t1){
double vx0=600*cos(M_PI/3.0);
double x = vx0*(1-exp(-gamma2*t1))/gamma2;
return(x);
}
//func¸~ao exata em y
double fy(double gamma3, double t2){
double vy0=600*sin(M_PI/3.0);
double y =(-10.0*t2/gamma3)+(1-exp(-gamma3*t2))*((gamma3*vy0+10.0)/(gamma3*gamma3));
return(y);
}
1.1 Variando h
Para analisar o impacto da mudanc¸a do passo h nos me´todos, e´ necessa´rio
criar um loop que ande com h e dentro deste loop os me´todos calculem a
soluc¸a˜o ate´ um ponto espec´ıfico para que o resultado possa ser comparado
com as soluc¸o˜es exatas. Apo´s a introduc¸a˜o deste loop, e definindo que os
me´todos calculara˜o ate´ o instante de tempo t = 50s, com γ = 0.05, o co´digo
fica conforme abaixo:
3
//NOME DO ARQUIVO: erros.c
//AUTOR: Fabio Rasera
#include <stdio.h>
#include <math.h>
double acc(double gamma1, double grav, double pos, double vel);
double fx(double gamma2, double t1);
double fy(double gamma3, double t2);
int main(void){
double t, E, E1, E0, theta, a0, a, h, x0, r[2], x, y, v0, y0, v[2], vx0, vy0, vx, vy, g[2];
double m1[2], m2[2], m3[2], m4[2], k1[2], k2[2], k3[2], k4[2], gamma;
//Dados iniciais
theta=(M_PI/3.0);
v0=600.0;
h=0.01;
gamma=0.05;
x0=0.0;
y0=0.0;
//Componentes em x e y da velocidade
vx0=v0*cos(theta);
vy0=v0*sin(theta);
//-----------RK2----------//
FILE *p;
p = fopen("RK2erro.txt","w");
for(h=pow(10,-6);h<=1;h*=10){
//Dados iniciais
g[0]=0.0; g[1]=10.0;
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
t=0;
do{
//Passo 1 // k1=f(t,x);
{
k1[0]= v[0];
k1[1]= v[1];
m1[0]= acc(gamma, g[0], r[0], k1[0]);
m1[1]= acc(gamma, g[1], r[1], k1[1]);
}
//Passo 2 // k2=f(t+h , x+k1*h*0.5);
{
k2[0]= v[0] + m1[0]*h*0.5;
k2[1]= v[1] + m1[1]*h*0.5;
m2[0]= acc(gamma, g[0], r[0] + k1[0]*h*0.5, k2[0]);
m2[1]= acc(gamma, g[1], r[1] + k1[1]*h*0.5, k2[1]);
}
r[0]= r[0] + k2[0]*h;
r[1]= r[1] + k2[1]*h;
v[0]= v[0] + m2[0]*h;
v[1]= v[1] + m2[1]*h;
t+=h;
}while(t<=50); //Repetir enquanto altura do proje´til for maior que zero
fprintf(p,"%.12e %.12e %.12e\n", h , ( fabs(r[0]-fx(gamma,t-h)) )/( fabs(fx(gamma,t-h)) ) , ( fabs(r[1]-fy(gamma,t-h)) )/( fabs(fy(gamma,t-h)) ) );
}
fclose(p);
//-----------RK4----------//
FILE *p1;
p1 = fopen("RK4erro.txt","w");
for(h=pow(10,-6);h<=1;h*=10){
//Reinicializac¸~ao dos dados
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
t=0;
4
do{
//Passo 1 // k1=f(t,x);
{
k1[0]= v[0];
k1[1]= v[1];
m1[0]= acc(gamma, g[0], r[0], k1[0]);
m1[1]= acc(gamma, g[1], r[1], k1[1]);
}
//Passo 2 // k2=f(t+0.5*h , x+0.5*k1*h);
{
k2[0]= v[0] + m1[0]*0.5*h;
k2[1]= v[1] + m1[1]*0.5*h;
m2[0]= acc(gamma, g[0], r[0]+0.5*k1[0]*h, k2[0]);
m2[1]= acc(gamma, g[1], r[1]+0.5*k1[1]*h, k2[1]);
}
//Passo 3 // k3=f(t+0.5*h, x+0.5*k2*h);
{
k3[0]= v[0] +m2[0]*0.5*h;
k3[1]= v[1] +m2[1]*0.5*h;
m3[0]= acc(gamma, g[0], r[0]+0.5*k2[0]*h, k3[0]);
m3[1]= acc(gamma, g[1], r[1]+0.5*k2[1]*h, k3[1]);
}
//Passo 4 // k4=f(t+h, x+k3*h);
{
k4[0]= v[0] +m3[0]*h;
k4[1]= v[1] +m3[1]*h;
m4[0]= acc(gamma, g[0], r[0]+k3[0]*h, k4[0]);
m4[1]= acc(gamma, g[1], r[1]+k3[1]*h, k4[1]);
}
//Completando o passo
{
r[0] = r[0] + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])*h/6.0;
r[1] = r[1] + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])*h/6.0;
v[0] = v[0] + (m1[0] + 2*m2[0] + 2*m3[0] + m4[0])*h/6.0;
v[1] = v[1] + (m1[1] + 2*m2[1] + 2*m3[1] + m4[1])*h/6.0;
t+=h;
}
}while(t<=50);
fprintf(p1,"%.12e%.12e %.12e\n", h , ( fabs(r[0]-fx(gamma,t-h)) )/( fabs(fx(gamma,t-h)) ) , ( fabs(r[1]-fy(gamma,t-h)) )/( fabs(fy(gamma,t-h)) ) );
}
fclose(p1);
//------Euler-Cromer------//
FILE *p2;
p2=fopen("Eulererro.txt","w");
for(h=pow(10,-6);h<=1;h*=10){
//Reinicializac¸~ao dos dados
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
t=0;
do{
v[0]+=acc(gamma, g[0], v[0], v[0])*h;
v[1]+=acc(gamma, g[1], v[1], v[1])*h;
r[0]+=v[0]*h;
r[1]+=v[1]*h;
t+=h;
}while(t<=50);
fprintf(p2,"%.12e %.12e %.12e\n", h , ( fabs(r[0]-fx(gamma,t-h)) )/( fabs(fx(gamma,t-h)) ) , ( fabs(r[1]-fy(gamma,t-h)) )/( fabs(fy(gamma,t-h)) ) );
}
fclose(p2);
5
return 0;
}
//func¸~ao acelerac¸~ao
double acc(double gamma1, double grav, double pos, double vel){
double a = -grav -gamma1*vel;
return(a);
}
//func¸~ao exata em x
double fx(double gamma2, double t1){
double vx0=600*cos(M_PI/3.0);
double x = vx0*(1-exp(-gamma2*t1))/gamma2;
return(x);
}
//func¸~ao exata em y
double fy(double gamma3, double t2){
double vy0=600*sin(M_PI/3.0);
double y =(-10.0*t2/gamma3)+(1-exp(-gamma3*t2))*((gamma3*vy0+10.0)/(gamma3*gamma3));
return(y);
}
Ao comparar os erros relativos de cada me´todo para as duas componen-
tes, nota-se que, com γ = 0.05, o erro dos me´todos RK2 e RK4 aumentam
na mesma proporc¸a˜o conforme o aumento de h, e o erro do Me´todo de Euler-
Cromer parece aumentar de forma semelhante, no entanto com um valor de
erro sempre maior que o dos outros me´todos. O erro na componente y se
demonstra um pouco maior que o erro na component x para um mesmo h.
O gra´fico a seguir mostra os resultados em escala logar´ıtmica.
 1e−09
 1e−08
 1e−07
 1e−06
 1e−05
 0.0001
 0.001
 0.01
 0.1
 1e−06 1e−05 0.0001 0.001 0.01 0.1 1
|x
n
u
m
−
x
e
x
a
to
| 
/ 
|x
e
x
a
to
|
h
Erro Relativo de Métodos Numéricos
Valores iniciais: γ=0.05, v0=600 m/s, θ0=pi/3, x0=0. y0=0
Euler−Cromer
RK2
RK4
Figura 1: Erro relativo na componente x da posic¸a˜o de um proje´til para diferentes
me´todos nume´ricos aplicados.
6
 1e−08
 1e−07
 1e−06
 1e−05
 0.0001
 0.001
 0.01
 0.1
 1
 1e−06 1e−05 0.0001 0.001 0.01 0.1 1
|y
n
u
m
−
y
e
x
a
to
| 
/ 
|y
e
x
a
to
|
h
Erro Relativo de Métodos Numéricos
Valores iniciais: γ=0.05, v0=600 m/s, θ0=pi/3, x0=0. y0=0
Euler−Cromer
RK2
RK4
Figura 2: Erro relativo na componente y da posic¸a˜o de um proje´til para diferentes
me´todos nume´ricos aplicados.
1.2 Variando γ
Para estudar o efeito do paraˆmetro γ, sera´ calculada a trajeto´ria do proje´til
para o me´todo RK4 conforme γ e´ variado. A escolha do incremento h
sera´ baseada no exerc´ıcio anterior. Pela Figura 1, veˆ-se que quanto menor
e´ h, menor o erro relativo em ambas as componentes. Pore´m, um fator
importante na escolha de h e´ que para h < 10−5 ja´ se notam efeitos na
compilac¸a˜o do programa, que se torna demorada ja´ que o tempo anda com
o incremento extremamente baixo. Assim, ja´ que com h = 10−3 o erro se
demonstra aceitavelmente baixo e na˜o aparenta acarretar em problemas de
compilac¸a˜o nem problemas de arredondamento, esse sera´ o h escolhido para
calcularmos os efeitos da constante de amortecimento. Os resultados para
as trajeto´rias do proje´til utilizando o me´todo RK4 e diferentes valores de γ
podem ser vistos no gra´fico abaixo.
7
 0
 2000
 4000
 6000
 8000
 10000
 12000
 14000
 0 5000 10000 15000 20000 25000 30000 35000
y
(m
)
x(m)
Posição de um projétil (RK4)
Valores iniciais: h=0.001, v0=600 m/s, θ0=pi/3, x0=0. y0=0
γ=0.0
γ=0.005
γ=0.01
γ=0.02
γ=0.04
γ=0.08
Figura 3: Posic¸a˜o de um proje´til conforme a constante de amortecimento e´ mo-
dificada no problema 1.
E´ evidente que, conforme se aumenta a constante de amortecimento,
menor e´ o caminho percorrido e a altura alcanc¸ada pelo proje´til.
2 Questa˜o 2
No pro´ximo problema a forc¸a de arrasto e´ modificada, resultando em um
par de equac¸o˜es diferenciais distinto do anterior:
x¨ = −γvx˙ (4)
y¨ = −mg − γvy˙ (5)
Abaixo o co´digo contendo os me´todos de Euler-Cromer, RK2 e RK4 para
a resoluc¸a˜o do novo problema:
//NOME DO ARQUIVO: Scherer.c
//AUTOR: Fabio Rasera
#include <stdio.h>
#include <math.h>
double acc(double gamma1, double grav, double pos, double velx, double vely, double velretorno);
int main(void){
double t, E, E1, E0, theta, a0, a, h, x0, r[2], x, y, v0, y0, v[2], vx0, vy0, vx, vy, g[2];
double m1[2], m2[2], m3[2], m4[2], k1[2], k2[2], k3[2], k4[2], gamma;
8
//Condic¸~oes iniciais
theta=(M_PI)/3.0;
h=0.01;
v0=600;
gamma=0.01;
x0=0.0;
y0=0.0;
vx0=v0*cos(theta);
vy0=v0*sin(theta);
g[0]=0.0; g[1]=10.0;
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
//------RK2------//
FILE *p;
p = fopen("RK2scherer.txt","w");
t=0;
do{
fprintf(p,"%.12e %.12e %.12e %.12e %.12e \n",r[0],r[1],v[0],v[1],t);
//Passo 1 // k1=f(t,x);
{
k1[0]=v[0];
k1[1]=v[1];
m1[0]= acc(gamma,g[0],r[0],k1[0],k1[1],k1[0]);
m1[1]= acc(gamma,g[1],r[1],k1[0],k1[1],k1[1]);
}
//Passo 2 // k2=f(t+h , x+k1*h*0.5);
{
k2[0]= v[0] + m1[0]*h*0.5;
k2[1]= v[1] + m1[1]*h*0.5;
m2[0]= acc(gamma, g[0], r[0] + k1[0]*h*0.5, k2[0], k2[1], k2[0]);
m2[1]= acc(gamma, g[1], r[1] + k1[1]*h*0.5, k2[0], k2[1], k2[1]);
}
r[0] = r[0] + k2[0]*h;
r[1] = r[1] + k2[1]*h;
v[0] = v[0] + m2[0]*h;
v[1] = v[1] + m2[1]*h;
t+=h;
}while(r[1]>=0);
fclose(p);
//------RK4------//
FILE *p1;
p1 = fopen("RK4scherer.txt","w");
//Reinicializando os dados
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
t=0;
do{
fprintf(p1,"%.12e %.12e %.12e %.12e %.12e \n",r[0],r[1],v[0],v[1],t);
//Passo 1 // k1=f(t,x);
{
k1[0]= v[0];
k1[1]= v[1];
m1[0]= acc(gamma, g[0], r[0], k1[0], k1[1], k1[0]);
m1[1]= acc(gamma, g[1], r[1], k1[0], k1[1], k1[1]);
}
//Passo 2 // k2=f(t+0.5*h , x+0.5*k1*h);
{
k2[0]= v[0] + m1[0]*0.5*h;
k2[1]= v[1] + m1[1]*0.5*h;
m2[0]= acc(gamma, g[0], r[0]+0.5*k1[0]*h, k2[0], k2[1], k2[0]);
9
m2[1]= acc(gamma, g[1], r[1]+0.5*k1[1]*h, k2[0], k2[1], k2[1]);
}
//Passo 3 // k3=f(t+0.5*h, x+0.5*k2*h);
{
k3[0]= v[0] +m2[0]*0.5*h;
k3[1]= v[1] +m2[1]*0.5*h;
m3[0]= acc(gamma, g[0], r[0]+0.5*k2[0]*h, k3[0], k3[1], k3[0]);
m3[1]= acc(gamma, g[1], r[1]+0.5*k2[1]*h, k3[0], k3[1], k3[1]);
}
//Passo 4 // k4=f(t+h, x+k3*h);
{
k4[0]= v[0] +m3[0]*h;
k4[1]= v[1] +m3[1]*h;
m4[0]= acc(gamma, g[0], r[0]+k3[0]*h, k4[0], k4[1], k4[0]);
m4[1]= acc(gamma, g[1], r[1]+k3[1]*h, k4[0], k4[1], k4[1]);
}
//Completando o passo
{
r[0] = r[0] + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])*h/6.0;
r[1] = r[1] + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])*h/6.0;
v[0] = v[0] + (m1[0] + 2*m2[0] + 2*m3[0] + m4[0])*h/6.0;
v[1] = v[1] + (m1[1] + 2*m2[1] + 2*m3[1] + m4[1])*h/6.0;
t+=h;
}
}while(r[1]>=0);
fclose(p1);
//----Euler-Cromer----//
FILE*p2;
p2=fopen("Eulerscherer.txt","w");
//Reinicializando os dados
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
t=0;
do{
fprintf(p2,"%.12e %.12e %.12e %.12e %.12e %.12e\n", r[0],r[1],v[0],v[1],h,t);
v[0]+=acc(gamma, g[0], r[0], v[0], v[1], v[0])*h;
v[1]+=acc(gamma, g[1], r[1], v[0], v[1], v[1])*h;
r[0]+=v[0]*h;
r[1]+=v[1]*h;
}while(r[1]>=0);
fclose(p2);
return 0;
}
//func¸~ao acelerac¸~ao
double acc(double gamma1,double grav, double pos, double velx, double vely, double velretorno){
double vel = sqrt(velx*velx + vely*vely);
double a = -grav -gamma1*vel*velretorno;
return(a);
}
2.1 Variando γ
Este problema foi resolvido utilizando o me´todo RK4. O gra´fico abaixo
demonstra o efeito da variac¸a˜o do coeficiente de amortecimento na soluc¸a˜o
do problema.
10
 0
 50
 100150
 200
 250
 300
 350
 400
 450
 0 50 100 150 200 250 300 350 400 450
y
 (
m
)
x (m)
Posição de um projétil (RK4)
Valores iniciais: h=0.001, v0=600 m/s, θ0=pi/3, x0=0. y0=0
γ=0.005
γ=0.01
γ=0.02
γ=0.04
γ=0.08
Figura 4: Posic¸a˜o de um proje´til conforme a constante de amortecimento e´ mo-
dificada no problema 2.
Assim como no problema anterior, nota-se a diminuic¸a˜o do alcance ho-
rizontal e vertical do proje´til conforme γ aumenta. Percebe-se, tambe´m,
que neste problema a forc¸a de arrasto atua com maior intensidade sobre o
proje´til, diminuindo ainda mais a distaˆncia e altura atingidas, mesmo que o
objeto tenha partido com velocidade inicial e aˆngulo iguais aos do problema
anterior.
3 Perguntas
3.1 Como avaliar o erro cometido quando na˜o conhecemos a
soluc¸a˜o exata?
Uma maneira de avaliar o erro dos me´todos de Runge Kutta e´ dividir o passo
h utilizado em 2 passos h/2, e definir que a diferenc¸a entre os resultados
obtido na˜o passe de um erro ma´ximo permitido. Assim, quando se testa
a diferenc¸a entre os resultados utilizando o passo h e os 2 passos h/2, o
programa pode modificar o passo anterior (h) para um passo novo que agora
passara´ a ser a metade, caso o erro viole o ma´ximo estipulado.
Ha´ outra maneira que exige menor capacidade de processamento do pro-
grama e funciona ta˜o bem quanto o me´todo anterior, que consiste basica-
mente em utilizar a precisa˜o de um meto´do de ordem m para estimar o erro
de um me´todo de ordem n < m.O problema evidente deste me´todo e´ que
11
se necessitarmos avaliar a precisa˜o do me´todo RK4, precisaremos de um
me´todo RK5. Este problema foi contornado por Fehlberg ao descobrir um
me´todo de quinta ordem que necessita do ca´lculo de 6 posic¸o˜es da func¸a˜o
dentro do intervalo h, mas que com uma combinac¸a˜o diferente de coeficien-
tes resulta em um me´todo de quarta ordem. Ou seja, a partir da escolha do
me´todo e coeficientes de Fehlberg, e´ poss´ıvel encontrar o erro corrente para
o me´todo RK4 ao compara´-lo com o RK5 e definir uma operac¸a˜o que decide
se o novo passo precisa ser mudado conforme o resultado deste erro.
Estes me´todos levam o nome de incrementos adaptativos, justamente
por adptarem o tamanho do incremento de acordo com a necessidade do
programa. Uma grande vantagem disso e´ que o incremento pode ser aumen-
tado quando a func¸a˜o e´ muito bem comportada e na˜o apresenta mudanc¸as
bruscas e diminu´ıdo caso ocorra o contra´rio. Assim, o incremento se torna
dinaˆmico e pode aumentar a velocidade de ca´lculo e precisa˜o do me´todo.
Utilizando o me´todo e coeficientes de Fehlberg para incrementos adap-
tativos, o novo valor do passo pode ser calculado da seguinte forma:
hnovo = h
∣∣∣∣ ��c
∣∣∣∣1/m (6)
Onde m e´ a ordem do me´todo de ordem mais alta, � e´ o erro ma´ximo
perimitido pelo usua´rio e �c e´ o erro corrente calculado pela diferenc¸a entre os
resultados dos me´todos de ordem m e n < m. Assim, quando o erro corrente
e´ maior que o erro ma´ximo permitido, h diminui, e se o erro corrente e´ menor,
h aumenta ate´ que o erro corrente chegue ao erro permitido, calculando a
soluc¸a˜o o mais rapidamente poss´ıvel dentro da precisa˜o desejada.
3.2 O que deveria ser modificado para que as coliso˜es com o
cha˜o fossem levadas em considerac¸a˜o?
Para que o programa continue a calcular a trajeto´ria do objeto levando
em considerac¸a˜o coliso˜es ela´sticas com o cha˜o e´ necessa´rio modificar a in-
formac¸a˜o da velocidade na componente y toda vez que a bola entra em
contato com o cha˜o. No entanto, e´ preciso tambe´m definir um limite para o
programa parar de calcular novas trajeto´rias, pois a velocidade inicial pode
chegar infinitamente pro´xima de zero. Minha escolha foi testar o desloca-
mento na componente x do objeto ate´ que ele seja muito pequeno em relac¸a˜o
a`s dimenso˜es do problema. No problema atual, quando esse deslocamento
for menor que 1mm podemos dizer que o objeto esta´ virtualmente parado.
Para demonstrar o resultado, sera´ utilizado o problema da Questa˜o 1, solu-
cionado com o me´todo RK4 e utilizando h = 0.01 e γ = 0.02. A seguir o
co´digo utilizado para a resoluc¸a˜o este problema.
12
//NOME DO ARQUIVO: colelas.c
//AUTOR: Fabio Rasera
#include <stdio.h>
#include <math.h>
double acc(double gamma1, double grav, double pos, double vel);
int main(void){
double xi, deltax, t, E, E1, E0, theta, a0, a, h, x0, r[2], x, y, v0, y0, v[2], vx0, vy0, vx, vy, g[2];
double m1[2], m2[2], m3[2], m4[2], k1[2], k2[2], k3[2], k4[2], gamma;
//Dados iniciais
theta=(M_PI/3.0);
v0=600.0;
h=0.01;
gamma=0.02;
x0=0.0; y0=0.0;
g[0]=0.0; g[1]=10.0;
//Componentes em x e y da velocidade
vx0=v0*cos(theta);
vy0=v0*sin(theta);
//-----------RK4----------//
FILE *p1;
p1 = fopen("RK4colelas.txt","w");
//Inicializac¸~ao dos dados
r[0]=x0; r[1]=y0;
v[0]=vx0; v[1]=vy0;
t=0;
xi=x0;
do{
fprintf(p1,"%.12e %.12e %.12e %.12e %.12e %.12e %.12e %.12e\n", r[0],r[1],v[0],v[1],h,t,fx(gamma,t),fy(gamma,t));
//Passo 1 // k1=f(t,x);
{
k1[0]= v[0];
k1[1]= v[1];
m1[0]= acc(gamma, g[0], r[0], k1[0]);
m1[1]= acc(gamma, g[1], r[1], k1[1]);
}
//Passo 2 // k2=f(t+0.5*h , x+0.5*k1*h);
{
k2[0]= v[0] + m1[0]*0.5*h;
k2[1]= v[1] + m1[1]*0.5*h;
m2[0]= acc(gamma, g[0], r[0]+0.5*k1[0]*h, k2[0]);
m2[1]= acc(gamma, g[1], r[1]+0.5*k1[1]*h, k2[1]);
}
//Passo 3 // k3=f(t+0.5*h, x+0.5*k2*h);
{
k3[0]= v[0] +m2[0]*0.5*h;
k3[1]= v[1] +m2[1]*0.5*h;
m3[0]= acc(gamma, g[0], r[0]+0.5*k2[0]*h, k3[0]);
m3[1]= acc(gamma, g[1], r[1]+0.5*k2[1]*h, k3[1]);
}
//Passo 4 // k4=f(t+h, x+k3*h);
{
k4[0]= v[0] +m3[0]*h;
k4[1]= v[1] +m3[1]*h;
m4[0]= acc(gamma, g[0], r[0]+k3[0]*h, k4[0]);
m4[1]= acc(gamma, g[1], r[1]+k3[1]*h, k4[1]);
}
//Completando o passo
{
r[0] = r[0] + (k1[0] + 2*k2[0] + 2*k3[0] + k4[0])*h/6.0;
r[1] = r[1] + (k1[1] + 2*k2[1] + 2*k3[1] + k4[1])*h/6.0;
v[0] = v[0] + (m1[0] + 2*m2[0] + 2*m3[0] + m4[0])*h/6.0;
v[1] = v[1] + (m1[1] + 2*m2[1] + 2*m3[1] + m4[1])*h/6.0;
t+=h;
}
13
if(r[1]<=0){
deltax = r[0]-xi;
xi=r[0];
v[1]=-v[1];
}else(deltax=1);
printf("%lf %lf %lf %lf\n",deltax, r[0], xi, r[1]);
}while(deltax>=pow(10,-3));
fclose(p1);
}
//func¸~ao acelerac¸~ao
double acc(double gamma1, double grav, double pos, double vel){
double a = -grav -gamma1*vel;
return(a);
}
Abaixo o gra´fico da trajeto´ria do objeto.
−1000
 0
 1000
 2000
 3000
 4000
 5000
 6000
 7000
 8000
 9000
 0 2000 4000 6000 8000 10000 12000 14000
y
(m
)
x(m)
Posição de um projétil com colisão elástica (RK4)
Valores iniciais: h=0.01, γ=0.02, v0=600 m/s, θ0=pi/3, x0=0. y0=0
Figura 5: Posic¸a˜o de um objeto lanc¸ado obliquamente com colisa˜o ela´stica ao
colidir com o cha˜o.
Percebe-se que a perda de energia em func¸a˜o da forc¸a resistiva e´ bastante
grande e que com γ = 0.02 o objeto alcanc¸a uma altura que e´ menor que a
metade da altura ma´xima inicial ja´ na segunda subida. No entanto, tambe´m
e´ nota´vel que o ganho em distaˆncia horizontal na˜o e´ trivial para este exemplo,
ja´ que o objeto percorre 25% a mais da distaˆncia em que acontece a primeira
colisa˜o.
14
	Exercício 1
	Variando h
	Variando 
	Questão 2
	Variando 
	Perguntas
	Como avaliar o erro cometido quando não conhecemos a solução exata?
	O que deveria ser modificado para que as colisões com o chão fossem levadas em consideração?