RK4 órbita de planeta, problemas no código

franciele.mendes
(usa Outra)

Enviado em 26/04/2015 - 14:22h

Tenho que criar um algoritmo para a órbita de um planeta, porém, sei que a força gravitacional se divide bidimensionalmente
Fx = (mMG / r³)*x
a = MGx / r³
Estou tendo problemas com esse r³ no programa, pois r = sqrt (x*x + y*y). Gostaria que alguém me ajudasse, se possível.
Também gostaria de saber se é necessário implementar o runge kutta com incrementos adaptativos, já que o erro vai ser muito grande em algumas partes da órbita (afélio e periélio) e como eu poderia fazer isso.
Abaixo segue o código:

#include <iostream>
#include <cmath>
#include <iomanip>

using namespace std;

int main()

{

float t;
float x = 0.0;
float y = 10.0;
float v_y = 0.0;
float v_x = 7.5;
float M = 2000.0;
float G = 1.0;
float tau = 0.001;
float t_final = 10.0;
float x1, v_x1, a_x1, x2, v_x2, a_x2, x3, v_x3, a_x3, x4, v_x4, a_x4;
float y1, v_y1, a_y1, y2, v_y2, a_y2, y3, v_y3, a_y3, y4, v_y4, a_y4;
float r, v;

while (t <= t_final)

{

t = t + tau;

x1 = x;
v_x1 = v_x;
a_x1 = (M*G*x) / (pow(x*x + y*y, 1.5));

x2 = x1 + v_x1*tau*0.5;
v_x2 = v_x1 + a_x1*tau*0.5;
a_x2 = (M*G*x1) / (pow(x1*x1 + y1*y1, 1.5));

x3 = x1 + v_x2*tau*0.5;
v_x3 = v_x1 + a_x2*tau*0.5;
a_x3 = (M*G*x2) / (pow(x2*x2 + y2*y2, 1.5));

x4 = x1 + v_x3*tau;
v_x4 = v_x1 + a_x3*tau;
a_x4 = (M*G*x4) / (pow(x4*x4 + y4*y4, 1.5));

y1 = y;
v_y1 = v_y;
a_y1 = (M*G*y) / (pow(x*x + y*y, 1.5));

y2 = y1 + v_y1*tau*0.5;
v_y2 = v_y1 + a_y1*tau*0.5;
a_y2 = (M*G*y1) / (pow(x1*x1 + y1*y1, 1.5));

y3 = y1 +v_y2*tau*0.5;
v_y3 = v_y1 + a_y2*tau*0.5;
a_y3 = (M*G*y2) / (pow(x2*x2 + y2*y2, 1,5));

y4 = y1 + v_y3*tau;
v_y4 = v_y1 + a_y3*tau;
a_y4 = (M*G*y4) / (pow(x4*x4 + y4*y4, 1.5));

x = x + ((v_x1 + 2*v_x2 + 2*v_x3 + v_x4)/6)*tau;

v_x = v_x + ((a_x1 + 2*a_x2 + 2*a_x3 + a_x4)/6)*tau;

y = y + ((v_y1 + 2*v_y2 + 2*v_y3 + v_y4)/6)*tau;

v_y = v_y + ((a_y1 + 2*a_y2 + 2*a_y3 + a_y4)/6)*tau;

r = sqrt (pow(x, 2) + pow(y, 2));
v = sqrt (pow(v_x, 2) + pow(v_y, 2));

//cout << x << "\t " << y << endl;
cout << r << "\t" << v << endl;
}


}
Obrigada.
Vou colocar aqui também o problema, caso alguém tenha dúvidas a respeito:
Imagine um sistema planetário composto por uma estrela massiva (mass M=2000) e um planeta (massa m=1). A única força que age no sistema é a gravidade (em módulo F = mMG/r^2, com G=1 neste caso e r sendo a distância entre o centro de massa dos corpos envolvidos) Da física básica, aprendemos que quando dois corpos interagem através de uma força central, o movimento deles se dá em um plano: portanto o problema pode ser tratado de forma bidimensional. Podemos, para facilitar os cálculos, colocar a estrela na origem do sistema (x=0, y=0) e, como a massa dela é muito maior do que a massa do planeta, podemos também desconsiderar o efeito da força de atração do planeta sobre a posição da estrela: portanto a posição da estrela pode ser, em uma primeira aproximação, considerada estática na origem.

O movimento do planeta pode ser separado em duas direções ortogonais, x e y. Assim, a força que atua sobre o planeta na direção x (Fx), por exemplo, pode ser descrita como:

Fx = ( mMG/r^2 ) * r_x = ( mMG/r^3 ) * x,

onde r_x é a projeção do vetor unitário da direção da força na direção x, que pode ser expandido como r_x = x / r, onde x é a projeção da distância entre o planeta e a estrela na direção x. O mesmo raciocínio vale para a direção y, e portanto:

Fy = ( mMG/r^2 ) * r_y = ( mMG/r^3 ) * y.

Assim, pede-se ao aluno que:

- Crie um algoritmo que calcula a posição (x e y) e velocidade (vx, vy) do planeta usando RK4 para as seguintes condições:

-> x0 = 0, y0 = 10, vy0=0, t0 = 0, dt=0.001, tf = 10

Para a velocidade inicial vx0, rode o programa com 3 valores diferentes: 7.5, sqrt(200) e 20.