Trabalho Órbitas com RK4 e Fehlberg (com programas)

Prévia do material em texto

Trabalho 3
Luan Bottin De Toni
IF-UFRGS
30 de junho de 2017
1 Introdução
Este trabalho consiste em quatro exercícios focados em problemas en-
volvendo forças gravitacionais. Primeiramente, será reproduzido o problema
encontrado no livro do Scherer [1] que envolve a órbita de um cometa, utili-
zando o método RK4 com passo temporal fixo e o método de Fehlberg com
incremento temporal adaptativo, ou seja, o programa irá ajustar o tamanho
do incremento temporal de acordo com o erro associado ao método naquele
passo específico.
Após analisada a diferença entre os métodos RK4 e Fehlberg, apenas
este último será usado para o restante dos exercícios. A próxima tarefa
será retratar a órbita de um dos planetas do sistema solar utilizando dados
reais disponíveis na literatura. Também é feito um exercício envolvendo um
sistema binário de estrelas, onde nenhuma delas se encontra na origem do
sistema, este caso pode ser feito quando as estrelas possuem massas iguais,
e quando possuem massas diferentes. Por último, será feito um problema de
três corpos, considerando o Sol na origem do sistema e calculando a órbita de
dois planetas, levando em consideração a aceleração provocada não apenas
pelo Sol, mas por todos corpos envolvidos.
2 Referencial Teórico
Alguns desses problemas possuem solução analítica, sendo possível com-
parar os resultados obtidos numericamente com os exatos.
A força gravitacional Newtoniana é uma força atrativa do tipo:
F (r) = −GMm
r2
(1)
onde M e m são as massas dos corpos considerados, G é a constante de
gravitação universal, e r é a distância entre os corpos. Como essa força é
conservativa há uma energia potencial associada a ela através da relação:
F (r) = −∇U(r) U(r) = −GMm
r
(2)
1
A energia total do sistema será a soma da energia potencial com a energia
cinética K = mv2/2:
E = K + U(r) (3)
Pela equação 1 e a segunda lei de Newton, sabemos que a aceleração do
corpo de massa m provocada pelo corpo de massa M será:
~a = −GM
r3
~r (4)
Sabe-se que o raio da órbita obedece à relação:
r(θ) =
l
�cos(θ) + 1
(5)
onde:
l =
L2
GMm
� =
√
1 +
2EL2
G2M2m
(6)
e L é o momentum angular, definido como:
L = mvr (7)
3 Exemplo do livro: órbita de um cometa
Buscou-se aqui reproduzir o exemplo encontrado no livro do Scherer [1].
Este programa calcula a órbita de um cometa com valores: GM = 1, r0 =
1000 e v0 = 0.02. A velocidade inicial se trata da velcidade tangencial do
cometa. Veremos no gráfico que o raio inicial é referente ao afélio, ou seja, a
maior distância em sua órbita, e a velocidade inicial é somente no eixo y.
Primeiramente o exemplo é reproduzido com método RK4 com passo
temporal h fixo, escolhido arbitrariamente como sendo h = 10. O programa
é rodado até completar uma órbita ou até o tempo de processamento chegar
ao limite estipulado no código. Calculou-se, a partir da equação 4, as velo-
cidades e posições nos eixos x e y para após transformá-las em coordenadas
polares r e θ. Também é plotado juntamente no gráfico a órbita exata do
cometa (calculada pela equação 5) a fim de assegurar o funcionamento do
método numérico.
O gráfico 1 mostra que as órbitas numérica e exata coincidem. Confir-
mando o que foi dito anteriormente, o afélio do cometa corresponde ao raio
inicial rmax = 1000, e seu periélio, calculado pelo programa, é rmin = 250.
A partir desses dados, e com o auxílio de algumas propriedades matemáticas
da elipse, podemos calcular a excentricidade da órbita � = 0, 6.
2
 500
 400
 300
 200
 100
 0
 100
 200
 300
 400
 500
 400 200 0 200 400 600 800 1000
 0 500 1000
y
x
RK4
RK4
Orbita exata
Figura 1: Órbita calculada pelo método RK4 com passo fixo.
Note que os pontos da órbita numérica se parecem com uma linha, de-
vido à sua proximidade. Isso significa que o programa está calculando muitos
pontos da órbita, provavelmente mais que o necessário. Além disso, há re-
giões onde não há a necessidade de calcular tantos pontos se comparado com
outras regiões da mesma órbita. O método de Fehlberg com passo temporal
variado nos deixa livres desta preocupação, pois este utiliza os métodos RK4
e RK5 para estimar o erro corrente (Ec = |xRK5 − xRK4|) em determinado
passo alterando o incremento temporal conforme a necessidade e baseando-se
em um valor pré-estabelecido de erro tolerável (neste caso, Et = 10
−8
); para
o caso de uma órbita, onde existe erro em x e y, será estabelecido como erro
corrente o maior entre esses dois.
A figura 2 mostra o gráfico da órbita com o método de Fehlberg. Nota-se
que as órbitas numérica e exata coincidem novamente, garantindo o funcio-
namento do método. Vemos, porém, que neste método o programa calcula
bem menos pontos se comparado ao RK4. Tendo em vista que foi estipulado
um erro tolerável consideravelmente pequeno para o programa, confirmamos,
assim, que o RK4 realmente fazia cálculos desnecessários para esse problema,
aumentando o tempo de processamento.
3
 500
 400
 300
 200
 100
 0
 100
 200
 300
 400
 500
 400 200 0 200 400 600 800 1000
 0 500 1000
y
x
Felhberg
Fehlberg
Orbita exata
Figura 2: Órbita calculada pelo método de Fehlberg.
Vamos agora analisar como o incremento temporal h se comporta no
programa; a suspeita inicial é que o erro corrente, e, por consequencia o
h, dependa do valor do raio da órbita, já que este depende está ligado à
velocidade. Assim, plotou-se um gráfico de r x h:
 0
 100
 200
 300
 400
 500
 600
 700
 800
 900
 200 300 400 500 600 700 800 900 1000
h
r
Felhberg − r x h
dados
Figura 3: Dependência do incremento temporal h com o raio r da órbita.
O gráfico 3 mostra claramente a dependência do h com o raio da órbita.
4
Quanto maior o raio, maior o incremento temporal usado, ou seja, o erro
corrente torna-se menor sendo possível usar um h maior nesses casos. Isso
quer dizer que o erro corrente é maior onde a velocidade do cometa é maior,
sendo necessário calcular mais ponto no periélio para o mesmo intervalo
de tempo no afélio. É possível notar no gráfico duas linhas, isso ocorre
porque o programa plotou os gráficos de uma órbita inteira, fazendo o cometa
percorrer o caminho entre seu raio máximo até o mínimo, e de volta do raio
mínimo até o máximo.
4 Órbita de um planeta ao redor do Sol
O objetivo desse exercício é tentar reproduzir a órbita de um dos planetas
do sistema solar, considerando o Sol na origem. Utilizaremos o método de
Fehlberg pelos motivos apresentados na seção anterior.
Esse tipo de problema lida com grandezas reais e, por ser de aspecto
astronômico, são de grande escala, tanto em massa, quanto em distância.
Porém, a constante G utilizada nos cálculos possui um valor bem pequeno
no S.I., sendo assim, procuramos redimensionar esta constante de acordo
com as unidades utilizadas nesse tipo de problema. Em astronomia, uma
unidade decorrente para a medida de distâncias é a unidade astrônomica
(UA) correspondente à distância da Terra até o Sol (149 597 870 700 metros).
A unidade de tempo será utilizada em anos; e a massa será expressa em
massas solares (Msol = 1, 989.10
30kg). Sendo assim:
G = 6.674287.10−11
m3
kg.s2
.
1UA3
(1, 496.1011m)3
.
1, 989.1030kg
1Msol
.
(31 557 600s)2
1ano2
G = 39, 488566
UA3
Msol ano2
(8)
Buscou-se reproduzir a órbita do planeta Saturno, o qual possui afélio
rmax = 9, 957 UA, velocidade mínima (no afélio) v0 ' 1, 92 UA/ano e massa
m ' 0, 000286 Msol. O erro tolerável utilizado pelo programa para calcular
o novo incremento temporal é de 10−10. Plotou-se a órbita do planeta com
essas condições iniciais (encontradas na referência [2]) juntamente com o
resultadoanalítico através da equação 5.
Vemos que o resultado numérico e analítico coincidem. Porém, se pesqui-
sarmos outros aspectos sobre o planeta podemos nos deparar com algumas
divergências como, por exemplo, o periélio, excentricidade e período do pla-
neta. O resultado numérico para o periélio é de aproximadamente 8, 6 UA
enquanto que o real encontrado na referência [2] é de 9, 195 UA. Da mesma
forma, a excentricidade da órbita mostra-se diferente, numericamente temos
o valor de 0, 0729 enquanto o valor real é de 0, 0565. Já o período é de
28, 14 anos para o método numérico, enquanto o real é de 29, 46 anos. Uma
5
 10
 8
 6
 4
 2
 0
 2
 4
 6
 8
 10
 10 8 6 4 2 0 2 4 6 8 10
 0 2 4 6 8 10y(U
A)
x(UA)
Orbita de Saturno
Numerica
Exata
Figura 4: Órbita de Saturno em torno do Sol.
explicação para essa discrepância nos resultados é de que estamos conside-
rando um sistema muito simples, com o Sol na origem do sistema, sendo este
o único a exercer força sobre o planeta, o que não condiz com a situação
física real do sistema solar.
-2e-03
-1e-03
-5e-04
0e+00
5e-04
1e-03
 0 5 10 15 20 25 30
E 
(U
A2
.
M
so
l/a
no
s2
)
t (anos)
Energia
Cinetica
Potencial
Total
Figura 5: Energia potencial, cinética e total de Saturno em sua órbita.
O gráfico 5 mostra como as energias potencial e cinética de Saturno
6
variam ao longo do tempo em sua órbita. O planeta possui energia cinética
máxima (e energia potencial mínima) quando sua velocidade é máxima, ou
seja, em seu periélio, na metade de seu período; da mesma forma, sua energia
cinética é mínima em seu afélio, no início e final da órbita. A energia total
constante está de acordo com a lei da conservação de energia; observando
seu valor negativo e lembrando das equações 2 e 3 conclui-se que a energia
potencial do planeta é maior que sua energia cinética, o que caracteriza um
estado ligado; se a energia cinética fosse maior que a potencial, o planeta
estaria acima de sua velocidade de escape, saindo de sua órbita.
5 Órbitas em um sistema binário de estrelas
Este exercício tem por objetivo reproduzir um sistema binário de estrelas
com massas iguais e com massas diferentes. O algoritmo deste programa é
semelhante ao anterior, calculando a órbita para cada estrela separadamente
considerando o centro de massa das duas estrelas:
CM =
(mara +mbrb)
ma +mb
(9)
ondema e ra são, respectivamente, a massa e a posição da estrela a; enquanto
mb e rb são as mesmas grandezas relacionadas à estrela b.
Sendo assim, o programa considera um corpo de massa M = ma +mb
localizado na origem e calcula as posições iniciais das estrelas a e b em relação
ao centro de massa para finalmente calcular suas órbitas pelo método de
Fehlberg, como nos exercícios anteriores. Vale notar que o erro corrente a
ser considerado neste caso será o maior entre as duas estrelas.
Para verificar o comportamente de duas estrelas de massas iguais utilizou-
se os seguintes parâmetros iniciais: ma = mb = 1 Msol; ra = −rb = 3 U.A.;
vya = −vyb = 2, 74 U.A./ano. Plotando um gráfico com tais parâmetros, e
um erro tolerável de 10−8, obtemos a figura 6.
Nota-se que as órbitas das duas estrelas são semelhantes, pois possuem
parâmetros iguais em módulo e, assim, estão sujeitas à mesma força em mó-
dulo. Apenas diferem na direção de suas órbitas; a estrela a, por possuir
velocidade inicial em y positiva, percorre a órbita no sentido anti-horário,
ao contrário da estrela b que possui velocidade tangencial inicial em y ne-
gativa e percorre a órbita em sentido horário. Seus períodos são iguais, de
aproximadamente 1, 6 anos.
7
 1.5
 1
 0.5
 0
 0.5
 1
 1.5
 3 2 1 0 1 2 3
 0 1 2 3y 
(U
A)
x (UA)
Sistema binario com massas iguais
Estrela a numerico
Estrela b numerico
Estrela a analitico
Estrela b analitico
Figura 6: Órbitas em um sistema binário com estrelas de massas iguais.
-2e+02
-2e+02
-1e+02
-5e+01
0e+00
5e+01
1e+02
2e+02
 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
E 
(U
A2
.
M
so
l/a
no
s2
)
t (anos)
Energia (massas iguais)
Cinetica (estrela a)
Potencial (estrela a)
Total (estrela a)
Cinetica (estrela b)
Potencial (estrela b)
Total (estrela b)
Figura 7: Energia das duas estrelas de massas iguais em um sistema binário.
A figura 7 mostra o gráfico das energias cinética, potencial e total das
estrelas em questão, sendo que as linhas referem-se à estrela a, enquanto os
pontos à estrela b. Nota-se que, assim como as órbitas, suas energias também
são semelhantes. A energia total constante nos diz que o sistema conserva
energia, como era esperado; enquanto seu valor negativo nos diz que se trata
de um sistema ligado, como explicado no exercício anterior.
Vamos agora considerar o caso em que as estrelas possuem massas dife-
rentes: ma = 1, 5 Msol; mb = 1 Msol. O restante dos parâmetros iniciais
8
permanecem os mesmos. Porém, note que o centro de massa não estará
equidistante das estrelas, diferentemente do problema anterior.
 1.5
 1
 0.5
 0
 0.5
 1
 1.5
 4 3 2 1 0 1 2 3
 0 1 2 3y 
(U
A)
x (UA)
Sistema binario com massas diferentes
Estrela a numerico
Estrela b numerico
Estrela a analitico
Estrela b analitico
Figura 8: Órbitas em um sistema binário com estrelas de massas diferentes.
O gráfico 8 mostra as órbitas das estrelas considerando o centro de massa
fixo na origem. A estrela a, por possuir maior massa iniciará a uma distância
menor do centro de massa (apesar de possuírem raios iniciais idênticos em
módulo) devido à equação 9; logo, essa estrela estará sujeita a uma intensi-
dade de força maior que a estrela b e, portanto, o módulo de sua velocidade
irá variar mais rapidamente, provocando uma órbita com período menor, já
que suas velocidades tangenciais iniciais também são iguais.
É interessante notar que a quantidade de pontos no gráfico para a estrela
b não é máxima em seu periélio como aconteceu no primeiro exercício deste
trabalho, e sim em regiões próximas de seu afélio. Isso é uma consequência
do modo como o erro corrente é considerado, sendo o maior entre os erros das
coordenadas das duas estrelas, ou seja, a região onde o programa utilizou um
incremento temporal pequeno pode ter sido por causa de um erro corrente
alto da trajetória da estrela a e refletiu-se nessa região específica da órbita
da outra estrela onde, a princípio, não haveria necessidade de um incremento
temporal tão pequeno.
A figura 9 mostra como as energias se comportam no caso de estrelas com
massas diferentes, sendo que as linhas representam as grandezas correspon-
dentes à estrela a e os pontos à estrela b. Vemos que nos dois casos a energia
total é constante, e esta mostra-se maior para a estrela b pois como as duas
estrelas iniciam com velocidades iguais, suas energias cinéticas também o
são, porém, sendo que a estrela b está mais distante do centro de massa, sua
energia potencial é maior. Note que, apesar de suas energias cinéticas iniciais
9
-8e+02
-6e+02
-4e+02
-2e+02
0e+00
2e+02
4e+02
6e+02
 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8
E 
(U
A2
.
M
so
l/a
no
s2
)
t (anos)
Energia (massas diferentes)
Cinetica (estrela a)
Potencial (estrela a)
Total (estrela a)
Cinetica (estrela b)
Potencial (estrela b)
Total (estrela b)
Figura 9: Energia das duas estrelas de massas diferentes em um sistema
binário.
serem iguais, tal energia atinge um valor muito maior para a estrela a, visto
que esta atinge uma velocidade maior em seu periélio por este ser menor
que o periélio da estrela b, adicionando o fato de que sua massa também é
maior. Fica evidente nesse gráfico a diferençanos períodos orbitais, sendo
que o período da estrela a é cerca de duas vezes maior que o da esrtela b.
6 Órbitas de dois planetas ao redor do Sol
Considerando o Sol na origem do sistema solar, será calculada a órbita de
dois planetas. Sendo que já foi calculada a órbita de Saturno neste trabalho,
será calculada novamente esta órbita a fim de comparar os resultados; o
outro planeta será Júpiter, por ser o planeta mais massivo do sistema solar
e, estando relativamente perto de Saturno, exerce maior força gravitacional
se comparado ao restante dos corpos do sistema solar.
Nessa situação, a aceleração sofrida pelos dois planetas não se origina
apenas do Sol, mas do outro planeta presente, ou seja, a aceleração resul-
tante no eixo x ou y será a soma da aceleração causada pelo Sol com a
do outro planeta. Para isso, o algoritmo foi escrito de modo que a função
"Fehlberg"receba, além dos parâmetros do planeta a ser calculado, também
os parâmetros do outro planeta para, assim, calcular suas distâncias rela-
tivas e a força exercida. A posição recebida pela função refere-se ao passo
anterior, ou seja, de t− h.
Novamente, será utilizado o método de Fehlberg com erro tolerável de
10−8. Os parâmetros de Saturno são os mesmos utilizados anteriormente, en-
quanto os de Júpiter, encontrados na mesma fonte, são: mj ' 0, 000954Msol;
raflio = 5, 369 UA; vafelio = 2, 62 UA/ano.
10
 10
 8
 6
 4
 2
 0
 2
 4
 6
 8
 10
 10 8 6 4 2 0 2 4 6 8 10
 0 2 4 6 8 10y 
(U
A)
x (UA)
Orbitas de Jupiter e Saturno
Jupiter
Saturno
Figura 10: Órbitas de Júpiter e Saturno em torno do Sol.
O gráfico 10 mostra as duas órbitas calculadas; o programa foi rodado
até que Saturno completasse sua órbita. Vemos que as atrações dos planetas
não é muito evidente, sendo que não provoca grande perturbação em suas
órbitas. A fim de verificar se, de fato, há alguma perturbação na órbita,
será comparada esta órbita de Saturno com a calculada no segundo exercício
deste trabalho para o mesmo planeta. Considerando que a órbita numérica
do segundo exercício coincindiu com a analítica, podemos usar esta última
para comparar com o presente resultado calculando a diferença relativa entre
os dois, ou seja, o módulo da diferença do raio da órbita de Saturno com
Júpiter presente e do raio da órbita analítica de Saturno dividido pelo módulo
deste último.
Fazendo o programa rodar até que t atinja 145 (anos), ou seja, Saturno
completa cerca de 5 órbitas; e aplicando escala logarítmica no eixo y, obteve-
se o gráfico 11 da diferença relativa das órbitas versus tempo. O gráfico
mostra que, de fato, há uma perturbação na órbita causada por Júpiter e em
alguns pontos a diferença dos raios pode ser pouco maior que 10%. Também
é interessante notar que a diferença em questão difere para cada uma das
órbitas de Saturno, isso acontece porque a órbita de Júpiter é menor e as
posições relativas dos dois planetas não é a mesma para todas órbitas. Por
exemplo, as condições iniciais colocam os dois planetas em seus afélios, no
eixo positivo de x, deixando-os muito próximo e fazendo com que a força exer-
cida por estes seja a maior possível, explicando por que o gráfico 11 começa
com um valor máximo para a diferença dos raios; quando Saturno completa
sua órbita, estando novamente em seu afélio, Júpiter já havia completado a
sua, estando em um ponto diferente do inicial.
11
 1e-06
 1e-05
 0.0001
 0.001
 0.01
 0.1
 1
 0 20 40 60 80 100 120 140
|S(
J) 
- S
|/|S
|
t (anos)
Diferenca orbita de Saturno
Figura 11: Diferença relativa entre o raio da órbita de Saturno com Júpiter
presente e do raio analítico da órbita de Saturno em função do tempo.
É razoável supor que, devido à atração gravitacional de Júpiter, a órbita
de Saturno irá diferir cada vez que este completar uma volta em torno do Sol.
Para ver isso acontecer, plotou-se um gráfico da órbita do planeta para um
tempo t = 3500 anos, aproximadamente 125 órbitas, com foco em seu afélio,
diferindo as primeiras órbitas das últimas pela cor da linha, como mostra o
gráfico 12.
 0.4
 0.2
 0
 0.2
 0.4
 9.86 9.88 9.9 9.92 9.94 9.96 9.98 10
 10y 
(U
A)
x (UA)
Orbita de Saturno em seu afelio
 0
 500
 1000
 1500
 2000
 2500
 3000
 3500
Figura 12: Órbita de Saturno em seu afélio ao longo de 3500 anos.
12
Vemos que há uma flutuação do ponto de afélio de cerca de 0, 06 UA, tal
diferença é provocada pela presença de Júpiter. E, se olharmos bem para as
cores das linhas, notaremos que esta diferença tende a se manter entre esses
pontos (aproximadamente entre 9, 9 e 9, 96) já que as cores mais amareladas,
que representam um instante de tempo bastante avançado, encontram-se es-
palhadas nessa faixa, indicando que a órbita não tende a sair desse intervalo.
7 Conclusão
Vimos no presente trabalho a vantagem de usar o método de Fehlberg,
comparando-o com o RK4. Tal método mostrou-se mais eficiente, no sentido
de fazer menor número de cálculos e, consequentemente, aprimorando o pro-
cessamento sem perda considerável de precisão, afinal, é possível estipular
um erro tolerável de modo que o programa adapte o incremento temporal
usado de acordo com este erro.
Logo, pôde-se utilizar tal método para resolver alguns problemas de me-
cânica gravitacional, reproduzindo a órbita de um planeta ao redor do Sol a
partir de dados reais; também mostrou-se como se comportam as órbitas de
duas estrelas em um sistema binário no caso em que suas massas são iguais
e no caso em que são diferentes. Por último, o problema de três corpos,
utilizando os planetas Júpiter e Saturno e comparando a órbita deste último
com a calculada anteriormente a fim de verificar o efeito que a presença de
Júpiter causa.
8 Códigos
8.1 Exercício 1 - Exemplo do Scherer com RK4
#include <stdio.h>
#include <math.h>
double fAcel(double rx, double ry);
void RK4(double h, double *px, double *py, double *pvx, double *pvy);
main () {
double x,y,vx,vy,r,theta;
double k,U;
double tmax=100000;
double t=0.0;
double GM=1;
double m=1.0;
13
double h=10.0;
int orb=0;
//condições iniciais
x=1000.0;
y=0.0;
vx=0.0;
vy=0.02;
FILE *p;
p=fopen("posicao.dat","w");
while ((t<=tmax) && (orb<2))
{
r=sqrt(x*x+y*y);
theta=atan2(y,x);
k = 0.5*m*(vx*vx + vy*vy);
U = -GM*m/r;
fprintf(p,"%.12e %.12e %.12e %.12e %.12e %.12e\n",theta,r,t,k,U,k+U);
//função RK4 recebe o end das variáveis e as atualiza diretamente
RK4(h,&x,&y,&vx,&vy);
//testa se completou a órbita
//a cada meia volta theta troca de sinal, 2 meias voltas = 1 volta -> orb=2
if (theta*atan2(y,x) < 0){
//printf("raio min/max: %lf\n",r);
orb++;
}
t+=h;
}
fclose(p);
}
void RK4 (double h, double *px, double *py, double *pvx, double *pvy)
{
14
double x=*px,y=*py,vx=*pvx,vy=*pvy;
double kx[5],ky[5],kvx[5],kvy[5];
double xm,ym,vmx,vmy;
//Bloco 1
kx[1] = h*vx;
kvx[1]= h*fAcel(x,y);
ky[1] = h*vy;
kvy[1]= h*fAcel(y,x);
//Bloco 2
//x
xm = x + 0.5*kx[1];
vmx= vx + 0.5*kvx[1];
kx[2] = h*vmx;
//y
ym = y + 0.5*ky[1];
vmy= vy + 0.5*kvy[1];
ky[2] = h*vmy;
//velocidades
kvx[2] = h*fAcel(xm,ym);
kvy[2] = h*fAcel(ym,xm);
//Bloco 3
//x
xm = x + 0.5*kx[2];
vmx = vx + 0.5*kvx[2];
kx[3] = h*vmx;
//y
ym = y + 0.5*ky[2];
vmy = vy + 0.5*kvy[2];
ky[3] = h*vmy;
// velocidades
kvx[3] = h*fAcel(xm,ym);
kvy[3] = h*fAcel(ym,xm);
//Bloco 4
//x
xm = x + kx[3];
vmx = vx + kvx[3];
kx[4] = h*vmx;
//y
ym = y + ky[3];
vmy = vy + kvy[3];
15
ky[4] = h*vmy;
// velocidades
kvx[4] = h*fAcel(xm,ym);
kvy[4] = h*fAcel(ym,xm);
//x e y
*px = x +(kx[1] + 2*kx[2] + 2*kx[3] + kx[4])/6.0;
*py = y + (ky[1] + 2*ky[2] + 2*ky[3] + ky[4])/6.0;
//vx e vy
*pvx = vx + (kvx[1] + 2*kvx[2] + 2*kvx[3] + kvx[4])/6.0;
*pvy = vy + (kvy[1] + 2*kvy[2] + 2*kvy[3] + kvy[4])/6.0;
}
double fAcel(double r1, double r2)
{
double GM=1.0;
double r=sqrt(r1*r1+r2*r2);
return (-GM*r1 / pow(r,3));
}
8.2 Exercício 1 - Exemplo do Scherer com método de Fehl-
berg
Utiliza a função "Fehlberg".
#include <stdio.h>
#include <math.h>
double _GM=1.0;
double fAcel(double r1, double r2);
#include "fehlberg.h"
main () {
double x,y,vx,vy,x0,y0,vx0,vy0,r,theta,xe,ye;
double k,U;
double tmax=100000;
double t=0.0;
double m=1.0;
double h=10,hnovo;
double S1=0.9, S2=2.0;
16
double Ecx,Ecy,Ec,Et;
int ok=1,orb=0;
//condições iniciais
x0=1000.0; x=x0;
y0=0.0; y=y0;
vx0=0.0; vx=vx0;
vy0=0.02; vy=vy0;
//erro tolerável
Et=pow(10,-8);
FILE *p;
p=fopen("posicaofel.dat","w");
while ((t<=tmax) && (orb<2))
{
//só irá imprimir se o erro corrente do passo anterior foi aceitável
if(ok==1){
r=sqrt(x0*x0 + y0*y0);
theta=atan2(y0,x0);
k = 0.5*m*(vx0*vx0 + vy0*vy0);
U = -_GM*m/r;
fprintf(p,"%.12e %.12e %.12e %.12e %.12e %.12e %.12e\n",theta,r,t,h,k,U,k+U);
}
//função Felhberg recebe o endereço das variáveis e altera seus valores diretamente. 'xe' e 'ye' serão utilizados para estimar o erro corrente em x e y
Fehlberg(h,&x,&y,&vx,&vy,&xe,&ye);
//estimando o erro corrente em x e y
Ecx = fabs(x - xe);
Ecy = fabs(y - ye);
//define o erro corrente como o maior entre o erro de x e de y
if (Ecx > Ecy) Ec=Ecx;
else Ec=Ecy;
if(Ec <= Et){ //erro aceitável
//atualiza as variáveis a serem utilizadas
x0 = x;
y0 = y;
17
vx0 = vx;
vy0 = vy;
//incrementa h
hnovo=h*S1*pow(Et/Ec, 0.2);
//Para evitar aumentos demasiados
if (hnovo > S2*h) hnovo=S2*h;
//testa se completou a órbita
//a cada meia volta theta troca de sinal, 2 meias voltas = 1 volta -> orb=2
if (theta*atan2(y0,x0) < 0) orb++;
ok=1;
t+=h;
}
else { //erro não aceitável
//repete os valores anteriores
x = x0;
y = y0;
vx = vx0;
vy = vy0;
//diminui h
hnovo=h*S1*pow(Et/Ec, 0.2);
ok=0;
}
h=hnovo;
}
fclose(p);
}
double fAcel(double r1, double r2)
{
double r;
r=sqrt(r1*r1+r2*r2);
return (-_GM*r1 / pow(r,3));
18
}
8.3 Função "Fehlberg
Função utilizada nos programas que utilizam o método de Fehlberg.
//para cálculo de uma órbita
//função felhberg utiliza a função fAcel. Recebe endereço das variáveis e as atualiza diretamente através de ponteiros
void Felhberg (double h, double *px, double *py, double *pvx, double *pvy, double *pxe, double *pye)
{
double x=*px,y=*py,vx=*pvx,vy=*pvy;
double kx[7],ky[7],kvx[7],kvy[7];
double xm,ym,vxm,vym;
//constantes utilizadas no método
double a2=1.0/5.0, a3=3.0/10.0, a4=3.0/5.0, a5=1.0, a6=7.0/8.0;
double b21=1.0/5.0, b31=3.0/40.0, b41=3.0/10.0, b51=-11.0/54.0, b61=1631.0/55296.0;
double b32=9.0/40.0, b42=-9.0/10.0, b52=5.0/2.0, b62=175.0/512.0;
double b43=6.0/5.0, b53=-70.0/27.0, b63=575.0/13824.0;
double b54=35.0/27.0, b64=44275.0/110592.0;
double b65=253.0/4096.0;
double c1=37.0/378.0, c2=0.0, c3=250.0/621.0, c4=125.0/594.0, c5=0.0, c6=512.0/1771.0;
double c1e=2825.0/27648.0, c2e=0.0, c3e=18575.0/48384.0, c4e=13525.0/55296.0, c5e=277.0/14336.0, c6e=1.0/4.0;
//Bloco 1
kx[1] = h*vx;
kvx[1]= h*fAcel(x,y);
ky[1] = h*vy;
kvy[1]= h*fAcel(y,x);
//Bloco 2
xm = x + b21*kx[1];
ym = y + b21*ky[1];
vxm = vx + b21*kvx[1];
vym = vy + b21*kvy[1];
kvx[2] = h*fAcel(xm , ym);
kvy[2] = h*fAcel(ym , xm);
19
kx[2] = h*vxm;
ky[2] = h*vym;
//Bloco 3
xm = x + b31*kx[1] + b32*kx[2];
ym = y + b31*ky[1] + b32*ky[2];
vxm = vx + b31*kvx[1] + b32*kvx[2];
vym = vy + b31*kvy[1] + b32*kvy[2];
kvx[3] = h*fAcel(xm, ym);
kvy[3] = h*fAcel(ym, xm);
kx[3] = h*vxm;
ky[3] = h*vym;
//Bloco 4
xm = x + b41*kx[1] + b42*kx[2] + b43*kx[3];
ym = y + b41*ky[1] + b42*ky[2] + b43*ky[3];
vxm = vx + b41*kvx[1] + b42*kvx[2] + b43*kvx[3];
vym = vy + b41*kvy[1] + b42*kvy[2] + b43*kvy[3];
kvx[4] = h*fAcel(xm, ym);
kvy[4] = h*fAcel(ym, xm);
kx[4] = h*vxm;
ky[4] = h*vym;
//Bloco 5
xm = x + b51*kx[1] + b52*kx[2] + b53*kx[3] + b54*kx[4];
ym = y + b51*ky[1] + b52*ky[2] + b53*ky[3] + b54*ky[4];
vxm = vx + b51*kvx[1] + b52*kvx[2] + b53*kvx[3] + b54*kvx[4];
vym = vy + b51*kvy[1] + b52*kvy[2] + b53*kvy[3] + b54*kvy[4];
kvx[5] = h*fAcel(xm, ym);
kvy[5] = h*fAcel(ym, xm);
kx[5] = h*vxm;
ky[5] = h*vym;
//Bloco 6
xm = x + b61*kx[1] + b62*kx[2] + b63*kx[3] + b64*kx[4] + b65*kx[5];
ym = y + b61*ky[1] + b62*ky[2] + b63*ky[3] + b64*ky[4] + b65*ky[5];
vxm = vx + b61*kvx[1] + b62*kvx[2] + b63*kvx[3] + b64*kvx[4] + b65*kvx[5];
20
vym = vy + b61*kvy[1] + b62*kvy[2] + b63*kvy[3] + b64*kvy[4] + b65*kvy[5];
kvx[6] = h*fAcel(xm, ym);
kvy[6] = h*fAcel(ym, xm);
kx[6] = h*vxm;
ky[6] = h*vym;
//x e y
*px = x + c1*kx[1] + c2*kx[2] + c3*kx[3] + c4*kx[4] + c5*kx[5] + c6*kx[6];
*py = y + c1*ky[1] + c2*ky[2] + c3*ky[3] + c4*ky[4] + c5*ky[5] + c6*ky[6];
//vx e vy
*pvx = vx + c1*kvx[1] + c2*kvx[2] + c3*kvx[3] + c4*kvx[4] + c5*kvx[5] + c6*kvx[6];
*pvy = vy + c1*kvy[1] + c2*kvy[2] + c3*kvy[3] + c4*kvy[4] + c5*kvy[5] + c6*kvy[6];
//cálculo de xe e ye para estimar o erro corrente
*pxe = x + c1e*kx[1] + c2e*kx[2] + c3e*kx[3] + c4e*kx[4] + c5e*kx[5] + c6e*kx[6];
*pye = y + c1e*ky[1] + c2e*ky[2] + c3e*ky[3] + c4e*ky[4] + c5e*ky[5] + c6e*ky[6];
}
8.4 Exercício 2 - Órbita de Saturno
Utiliza a função "Fehlberg"
#include <stdio.h>
#include <math.h>
double _GM=39.488566; //UA/a
double fAcel(double r1, double r2);
#include "fehlberg.h"
main () {
double x,y,vx,vy,x0,y0,vx0,vy0,r,theta,xe,ye;
double k,U;
double tmax=50; //anos
double t=0.0;
double m=0.000285741579; //em massas solares
double h=0.02,hnovo; //em anos
double S1=0.9, S2=2.0;
double Ecx,Ecy,Ec,Et;
int ok=1, orb=0;
21
//condições iniciais para Saturno
x0=9.957; x=x0; //em UA
y0=0.0; y=y0;
vx0=0.0; vx=vx0;
vy0=1.91750390; vy=vy0; //em UA/ano
//erro tolerável
Et=pow(10,-10);
FILE *p;
p=fopen("posicao.dat","w");
while ((t<=tmax) && (orb<2))
{
//só irá imprimir se o erro corrente do passo anterior foi aceitável
if(ok==1){
r=sqrt(x0*x0 + y0*y0);
theta=atan2(y0,x0);
k = 0.5*m*(vx0*vx0 + vy0*vy0);
U = -_GM*m/r;
fprintf(p,"%.12e %.12e %.12e %.12e %.12e %.12e %.12e\n",theta,r,t,h,k,U,k+U);
}
//função Fehlberg recebe o endereço das variáveis e altera seus valores diretamente. 'xe' e 'ye' serão utilizados para estimar o erro corrente em x e y
Fehlberg(h,&x,&y,&vx,&vy,&xe,&ye);
//estimando o erro corrente em x e y
Ecx = fabs(x - xe);
Ecy = fabs(y - ye);
//define o erro corrente como o maior entre o erro de x e de y
if (Ecx > Ecy) Ec=Ecx;
else Ec=Ecy;
if(Ec <= Et){ //erro aceitável
//atualiza as variáveis a serem utilizadas
x0 = x;
y0 = y;
vx0 = vx;
vy0 = vy;
22
//testa se completou a órbita
//a cada meia volta theta troca de sinal, 2 meias voltas = 1 volta -> orb=2
if (theta*atan2(y0,x0) < 0) orb++;
//incrementa h
hnovo=h*S1*pow(Et/Ec, 0.2);
//Para evitar aumentos demasiados
if (hnovo > S2*h) hnovo=S2*h;
ok=1;
t+=h;
}
else { //erro não aceitável
//repete os valores anteriores
x = x0;
y = y0;
vx = vx0;
vy = vy0;
//diminui h
hnovo=h*S1*pow(Et/Ec, 0.2);
ok=0;
}
h=hnovo;
}
fclose(p);
}
double fAcel(double r1, double r2)
{
double r;
r=sqrt(r1*r1+r2*r2);
return (-_GM*r1 / pow(r,3));
}
23
8.5 Exercício 3 - Sistema binário
Utiliza a função "Fehlberg"
#include <stdio.h>
#include <math.h>
double _G=39.488566;//UA/M_sol.a
double _M;
double fAcel(double r1, double r2);
#include "fehlberg.h"
main () {
double xa,xb,ya,yb,vxa,vxb,vya,vyb;
double x0a,y0a,vx0a,vy0a,x0b,y0b,vx0b,vy0b;
double r_a,theta_a,xea,yea;
double r_b,theta_b,xeb,yeb;
double ma,mb,CM;
double k,U;
double tmax=1.9; //anos
double t=0.0;
double h=0.02,hnovo; //em anos
double S1=0.9, S2=2.0;
double Ecx,Ecy,Ec,Eca,Ecb,Et;
int ok=1;
//condições iniciais do corpo a
ma=1.5; //em massas solares
x0a=3.0; xa=x0a; //em UA
y0a=0.0; ya=y0a;
vx0a=0.0; vxa=vx0a;
vy0a=2.74; vya=vy0a; //em UA/ano
//condições iniciais do corpo b
mb=1.0; //em massas solares
x0b=-3.0; xb=x0b; //em UA
y0b=0.0; yb=y0b;
vx0b=0.0; vxb=vx0b;
vy0b=-2.74; vyb=vy0b; //em UA/ano
//erro tolerável
24
Et=pow(10,-8);
//calcula centro de massa e define as posições dos corpos a e b de modo que o CM fique na origem
CM = (ma*xa + mb*xb)/(ma+mb);
xa = xa-CM; x0a=xa;
xb = xb-CM; x0b=xb;
_M = ma+mb;
FILE *pa,*pb;
pa=fopen("posicao_a2.dat","w");
pb=fopen("posicao_b2.dat","w");
while (t<=tmax)
{
//só irá imprimir se o erro corrente do passo anterior foi aceitável
if(ok==1){
r_a=sqrt(x0a*x0a + y0a*y0a);
theta_a=atan2(y0a,x0a);
k = 0.5*ma*(vx0a*vx0a + vy0a*vy0a);
U = -_G*_M*ma/r_a;
fprintf(pa,"%.12e %.12e %.12e %.12e %.12e %.12e %.12e\n",theta_a,r_a,t,h,k,U,k+U);
r_b=sqrt(x0b*x0b + y0b*y0b);
theta_b=atan2(y0b,x0b);
k = 0.5*mb*(vx0b*vx0b + vy0b*vy0b);
U = -_G*_M*mb/r_b;
fprintf(pb,"%.12e %.12e %.12e %.12e %.12e %.12e %.12e\n",theta_b,r_b,t,h,k,U,k+U);
}
//função Fehlberg recebe o endereço das variáveis e altera seus valores diretamente. 'xe' e 'ye' serão utilizados para estimar o erro corrente em x e y
Fehlberg(h,&xa,&ya,&vxa,&vya,&xea,&yea);
Fehlberg(h,&xb,&yb,&vxb,&vyb,&xeb,&yeb);
//Erro corpo a
//estimando o erro corrente em x e y
Ecx = fabs(xa - xea);
Ecy = fabs(ya - yea);
//define o erro corrente em a como o maior entre o erro de x e de y
25
if (Ecx > Ecy) Eca=Ecx;
else Eca=Ecy;
//Erro corpo b
//estimando o erro corrente em x e y
Ecx = fabs(xb - xeb);
Ecy = fabs(yb - yeb);
//define o erro corrente em b como o maior entre o erro de x e de y
if (Ecx > Ecy) Ecb=Ecx;
else Ecb=Ecy;
//define o erro corrente como o maior entre a e b
if (Eca > Ecb) Ec=Eca;
else Ec=Ecb;
if(Ec <= Et){ //erro aceitável
//atualiza as variáveis a serem utilizadas
x0a = xa; x0b = xb;
y0a = ya; y0b = yb;
vx0a = vxa; vx0b = vxb;
vy0a = vya; vy0b = vyb;
//incrementa h
hnovo=h*S1*pow(Et/Ec, 0.2);
//Para evitar aumentos demasiados
if (hnovo > S2*h) hnovo=S2*h;
ok=1;
t+=h;
}
else { //erro não aceitável
//repete os valores anteriores
xa = x0a; xb = x0b;
ya = y0a; yb = y0b;
vxa = vx0a; vxb = vx0b;
vya = vy0a; vyb = vy0b;
//diminui h
hnovo=h*S1*pow(Et/Ec, 0.2);
ok=0;
}
26
h=hnovo;
}
fclose(pa);
fclose(pb);
}
double fAcel(double r1, double r2)
{
double r;
r=sqrt(r1*r1+r2*r2);
return (-_G*_M*r1 / pow(r,3));
}
8.6 Exercício 4 - Sistema de 3 corpos
Utiliza a função "Fehlberg2"
#include <stdio.h>
#include <math.h>
double _G=39.488566; //UA/M_sol.a
double fAcel(double r1, double r2, double M);
double fSatAnalitica(double t);
#include "fehlberg2.h"
main () {
double xa,xb,ya,yb,vxa,vxb,vya,vyb;
double x0a,y0a,vx0a,vy0a,x0b,y0b,vx0b,vy0b;
double r_a,theta_a,xea,yea;
double r_b,theta_b,xeb,yeb;
double ma,mb;
double tmax=28.5; //anos
double t=0.0;
double h=0.02,hnovo; //em anos
double S1=0.9, S2=2.0;
double Ecx,Ecy,Ec,Eca,Ecb,Et,erro;
27
int ok=1;
//condições iniciais de Saturno
ma=0.000285741579; //em massas solares
x0a=9.957; xa=x0a; //em UA
y0a=0.0; ya=y0a;
vx0a=0.0; vxa=vx0a;
vy0a=1.9175039; vya=vy0a; //em UA/ano
//condições iniciais de Jupiter
mb=0.000954382278; //em massas solares
x0b=5.369; xb=x0b; //em UA
y0b=0.0; yb=y0b;
vx0b=0.0; vxb=vx0b;
vy0b=2.62421199; vyb=vy0b; //em UA/ano
//erro tolerável
Et=pow(10,-8);
FILE *pa,*pb;
pa=fopen("posicao_S.dat","w");
pb=fopen("posicao_J.dat","w");
while (t<=tmax)
{
//só irá imprimir se o erro corrente do passo anterior foi aceitável
if(ok==1){
//Saturno
r_a=sqrt(x0a*x0a + y0a*y0a);
theta_a=atan2(y0a,x0a);
erro=fabs(r_a-fSatAnalitica(t))/fabs(fSatAnalitica(t));
fprintf(pa,"%.12e %.12e %.12e %.12e %.12e\n",theta_a,r_a,t,h,erro);
//Jupiter
r_b=sqrt(x0b*x0b + y0b*y0b);
theta_b=atan2(y0b,x0b);
fprintf(pb,"%.12e %.12e %.12e %.12e\n",theta_b,r_b,t,h);
}
//função Fehlberg recebe o endereço das variáveis e altera seus valores diretamente. 'xe' e 'ye' serão utilizados para estimar o erro corrente em x e y
//também recebe os parâmetros do outro planeta em questão para calcular a aceleração referente a ele
28
Fehlberg2(h,&xa,&ya,&vxa,&vya,x0b,y0b,mb,&xea,&yea); //saturno
Fehlberg2(h,&xb,&yb,&vxb,&vyb,x0a,y0a,ma,&xeb,&yeb); //jupiter
//Erro corpo a
//estimando o erro corrente em x e y
Ecx = fabs(xa - xea);
Ecy = fabs(ya - yea);
//define o erro corrente em a como o maior entre o erro de x e de y
if (Ecx > Ecy) Eca=Ecx;
else Eca=Ecy;
//Erro corpo b
//estimando o erro corrente em x e y
Ecx = fabs(xb - xeb);
Ecy = fabs(yb - yeb);
//define o erro corrente em b como o maior entre o erro de x e de y
if (Ecx > Ecy) Ecb=Ecx;
else Ecb=Ecy;
//define o erro corrente como o maior entre a e b
if (Eca > Ecb) Ec=Eca;
else Ec=Ecb;
if(Ec <= Et){ //erro aceitável
//atualiza as variáveis a serem utilizadas
x0a = xa; x0b = xb;
y0a = ya; y0b = yb;
vx0a = vxa; vx0b = vxb;
vy0a = vya; vy0b = vyb;
//incrementa h
hnovo=h*S1*pow(Et/Ec, 0.2);
//Para evitar aumentos demasiados
if (hnovo > S2*h) hnovo=S2*h;
ok=1;
t+=h;
}
else { //erro não aceitável
//repete os valores anteriores
xa = x0a; xb = x0b;
ya = y0a; yb = y0b;
vxa = vx0a; vxb = vx0b;
29
vya = vy0a; vyb = vy0b;
//diminui h
hnovo=h*S1*pow(Et/Ec, 0.2);
ok=0;
}
h=hnovo;
}
fclose(pa);
fclose(pb);
}
double fAcel(double r1, double r2, double M)
{
double r;
r=sqrt(r1*r1+r2*r2);
return (-_G*M*r1 / pow(r,3));
}
double fSatAnalitica(double t)
{
double x0=9.957;
double v0=1.9175039;
double L=x0*v0;
double E=v0*v0/2 - _G/x0;
double ec=sqrt(1 + 2*E*L*L/(_G*_G));
double a=_G/(2*E);
double l=a*(1-ec*ec);
return (-l/(1+ec*cos(t)));
// return( (x0/(v0*v0*x0-2))*(-pow(v0,4)*x0*x0+2*x0*v0*v0)*(1/(1+cos(t)*sqrt(1+x0*x0*pow(v0,4)-2*x0*v0*v0))));
}
30
8.7 Função Fehlberg2
Função utilizada no problema de 3 corpos.
//para cálculo de uma órbita com dois corpos
//calcula a aceleração causada pelo Sol e pelo outro planeta, recebe as coordenadas referentes ao outro planeta juntamente com sua massa
//Utiliza a função fAcel enviando como parâmetro a massa e as coordenadas relativas ao outro corpo (sol ou planeta)
//fAcel(x,y,1) é a aceleração referente ao sol
//fAcel(x-X,y-Y,m) é a aceleração referente ao outro planeta
void Fehlberg2 (double h, double *px, double *py, double *pvx, double *pvy, double X, double Y, double m, double *pxe, double *pye)
{
double x=*px,y=*py,vx=*pvx,vy=*pvy;
double kx[7],ky[7],kvx[7],kvy[7];
double xm,ym,vxm,vym;
//constantes utilizadas no método
double a2=1.0/5.0, a3=3.0/10.0, a4=3.0/5.0, a5=1.0, a6=7.0/8.0;
double b21=1.0/5.0, b31=3.0/40.0, b41=3.0/10.0, b51=-11.0/54.0, b61=1631.0/55296.0;
double b32=9.0/40.0, b42=-9.0/10.0, b52=5.0/2.0, b62=175.0/512.0;
double b43=6.0/5.0, b53=-70.0/27.0, b63=575.0/13824.0;
double b54=35.0/27.0, b64=44275.0/110592.0;
double b65=253.0/4096.0;
double c1=37.0/378.0, c2=0.0, c3=250.0/621.0, c4=125.0/594.0, c5=0.0, c6=512.0/1771.0;
double c1e=2825.0/27648.0, c2e=0.0, c3e=18575.0/48384.0,c4e=13525.0/55296.0, c5e=277.0/14336.0, c6e=1.0/4.0;
//Bloco 1
kx[1] = h*vx;
kvx[1]= h*(fAcel(x,y,1.0) + fAcel(x-X,y-Y,m));
ky[1] = h*vy;
kvy[1]= h*(fAcel(y,x,1.0) + fAcel(y-Y,x-X,m));
//Bloco 2
xm = x + b21*kx[1];
ym = y + b21*ky[1];
vxm = vx + b21*kvx[1];
vym = vy + b21*kvy[1];
kvx[2] = h*(fAcel(xm,ym,1.0) + fAcel(xm-X,ym-Y,m));
kvy[2] = h*(fAcel(ym,xm,1.0) + fAcel(ym-Y,xm-X,m));
kx[2] = h*vxm;
ky[2] = h*vym;
31
//Bloco 3
xm = x + b31*kx[1] + b32*kx[2];
ym = y + b31*ky[1] + b32*ky[2];
vxm = vx + b31*kvx[1] + b32*kvx[2];
vym = vy + b31*kvy[1] + b32*kvy[2];
kvx[3] = h*(fAcel(xm,ym,1.0) + fAcel(xm-X,ym-Y,m));
kvy[3] = h*(fAcel(ym,xm,1.0) + fAcel(ym-Y,xm-X,m));
kx[3] = h*vxm;
ky[3] = h*vym;
//Bloco 4
xm = x + b41*kx[1] + b42*kx[2] + b43*kx[3];
ym = y + b41*ky[1] + b42*ky[2] + b43*ky[3];
vxm = vx + b41*kvx[1] + b42*kvx[2] + b43*kvx[3];
vym = vy + b41*kvy[1] + b42*kvy[2] + b43*kvy[3];
kvx[4] = h*(fAcel(xm,ym,1.0) + fAcel(xm-X,ym-Y,m));
kvy[4] = h*(fAcel(ym,xm,1.0) + fAcel(ym-Y,xm-X,m));
kx[4] = h*vxm;
ky[4] = h*vym;
//Bloco 5
xm = x + b51*kx[1] + b52*kx[2] + b53*kx[3] + b54*kx[4];
ym = y + b51*ky[1] + b52*ky[2] + b53*ky[3] + b54*ky[4];
vxm = vx + b51*kvx[1] + b52*kvx[2] + b53*kvx[3] + b54*kvx[4];
vym = vy + b51*kvy[1] + b52*kvy[2] + b53*kvy[3] + b54*kvy[4];
kvx[5] = h*(fAcel(xm,ym,1.0) + fAcel(xm-X,ym-Y,m));
kvy[5] = h*(fAcel(ym,xm,1.0) + fAcel(ym-Y,xm-X,m));
kx[5] = h*vxm;
ky[5] = h*vym;
//Bloco 6
xm = x + b61*kx[1] + b62*kx[2] + b63*kx[3] + b64*kx[4] + b65*kx[5];
ym = y + b61*ky[1] + b62*ky[2] + b63*ky[3] + b64*ky[4] + b65*ky[5];
vxm = vx + b61*kvx[1] + b62*kvx[2] + b63*kvx[3] + b64*kvx[4] + b65*kvx[5];
vym = vy + b61*kvy[1] + b62*kvy[2] + b63*kvy[3] + b64*kvy[4] + b65*kvy[5];
kvx[6] = h*(fAcel(xm,ym,1.0) + fAcel(xm-X,ym-Y,m));
32
kvy[6] = h*(fAcel(ym,xm,1.0) + fAcel(ym-Y,xm-X,m));
kx[6] = h*vxm;
ky[6] = h*vym;
//x e y
*px = x + c1*kx[1] + c2*kx[2] + c3*kx[3] + c4*kx[4] + c5*kx[5] + c6*kx[6];
*py = y + c1*ky[1] + c2*ky[2] + c3*ky[3] + c4*ky[4] + c5*ky[5] + c6*ky[6];
//vx e vy
*pvx = vx + c1*kvx[1] + c2*kvx[2] + c3*kvx[3] + c4*kvx[4] + c5*kvx[5] + c6*kvx[6];
*pvy = vy + c1*kvy[1] + c2*kvy[2] + c3*kvy[3] + c4*kvy[4] + c5*kvy[5] + c6*kvy[6];
//cálculo de xe e ye para estimar o erro corrente
*pxe = x + c1e*kx[1] + c2e*kx[2] + c3e*kx[3] + c4e*kx[4] + c5e*kx[5] + c6e*kx[6];
*pye = y + c1e*ky[1] + c2e*ky[2] + c3e*ky[3] + c4e*ky[4] + c5e*ky[5] + c6e*ky[6];
}
Referências
[1] CLAUDIO SCHERER Métodos Computacionais da Física, Editora Li-
vraria da Física, 2
a
ed., São Paulo, 2010.
[2] NASA Planetary Fact Sheet - Metric. Disponível em: http:
//nssdc.gsfc.nasa.gov/planetary/factsheet/index.html. Aces-
sado em 29/04/2016.
33
	Introdução
	Referencial Teórico
	Exemplo do livro: órbita de um cometa
	Órbita de um planeta ao redor do Sol
	Órbitas em um sistema binário de estrelas
	Órbitas de dois planetas ao redor do Sol
	Conclusão
	Códigos
	Exercício 1 - Exemplo do Scherer com RK4
	Exercício 1 - Exemplo do Scherer com método de Fehlberg
	Função "Fehlberg
	Exercício 2 - Órbita de Saturno
	Exercício 3 - Sistema binário
	Exercício 4 - Sistema de 3 corpos
	Função Fehlberg2