## Two Body problem with RK4 in C++

I want to output an excel file with the results of the trajectories of a two body problem, with initial position and velocity. But my program is not compiling. Any suggestions/problems that you can see?

Code:
#include <iostream>
#include <vector>
#include <fstream>

using namespace std;

void RK4_func(double,double,double,double,double,double,double,double,double,double,double,double);

int main(){

//initialise variables
double time_end,time_step,mass1,x1x,y1y,v1x,v1y,mass2,x2x,y2y,v2x,v2y;

//Get the value
cout << "Enter Simulation Duration in seconds" << endl;
cin >> time_step;
cout << "\n";

cout << "Enter Precision (Time step)" << endl;
cin >> time_end;
cout << "\n";

cout << "Enter mass of body 1 with 10^20kg scale" << endl;
cin >> mass1;
cout << "\n";

cout << "Enter initial body 1 x position" << endl;
cin >> x1x;
cout << "\n";

cout << "Enter initial body 1 y position" << endl;
cin >> y1y;
cout << "\n";

cout << "Enter initial body 1 x velocity" << endl;
cin >> v1x;
cout << "\n";

cout << "Enter initial body 1 x velocity" << endl;
cin >> v1y;
cout << "\n";

cout << "Enter mass of body 2 with 10^20kg scale" << endl;
cin >> mass2;
cout << "\n";

cout << "Enter initial body 2 x position" << endl;
cin >> x2x;
cout << "\n";

cout << "Enter initial body 2 y position" << endl;
cin >> y2y;
cout << "\n";

cout << "Enter initial body 2 x velocity" << endl;
cin >> v2x;
cout << "\n";

cout << "Enter initial body 2 y velocity" << endl;
cin >> v2y;
cout << "\n";

//Pass values int RK4_func
mass1=mass1*10^20
mass2=mass2*10^20
RK4_func(time_end,time_step;mass1;x1x,y1y,v1x,v1y,mass2,x2x,y2y,v2x,v2y);

//Output data

return 0;
}

//Assume function  enters data in SI units

void RK4_func (double timeend, double steps, double m1, double x1, double y1, double vx1, double vy1, double m2, double x2, double y2, double vx2,double vy2)
{

//Excel file creation

ofstream outfile1 ("RK4.xls");
outfile1<<"time"<<'\t'<<"x1"<<'\t'<<"y1"<<'\t'<<"vx1"<<'\t'<<"vy1"<<'\t'<<"x2"<<'\t'<<"y2"<<'\t'<<"vx2"<<'\t'<<"vy2"<<endl;
outfile1<<0<<'\t'<<x1<<'\t'<<y1<<'\t'<<vx1<<'\t'<<vy1<<'\t'<<x2<<'\t'<<y2<<'\t'<<vx2<<'\t'<<vy2<<endl;

constant double G = 6.67E-11;
double datalength = 1+(timeend/steps);
double data [8][datalength];

double t=0;
data[0][0]=x1;
data[1][0]=y1;
data[2][0]=x2;
data[3][0]=y2;
data[4][0]=vx1;
data[5][0]=vy1;
data[6][0]=vx2;
data[7][0]=vy2;

for (int i=0; i<datalength; i++)
{
vector<double> interdata(12), k1 (8), k2 (8), k3 (8), k4 (8);

// let 0~3 indicate position x1y1x2y2 4~7 indicate vx1vy1vx2vy2 8~11 indicate ax1ay1ax2ay2

//acceleration x1 y1 x2 y2

interdata[0]=data[0][i];
interdata[1]=data[1][i];
interdata[2]=data[2][i];
interdata[3]=data[3][i];
interdata[4]=data[4][i];
interdata[5]=data[5][i];
interdata[6]=data[6][i];
interdata[7]=data[7][i];
interdata[8] = G*m2*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[2][i]-data[0][i]);
interdata[9] = G*m2*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[3][i]-data[1][i]);
interdata[10] = G*m1*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[0][i]-data[2][i]);
interdata[11] = G*m1*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[1][i]-data[3][i]);

//k1 is Euler Step

k1[0] = steps*interdata[8];
k1[1] = steps*interdata[9];
k1[2] = steps*interdata[10];
k1[3] = steps*interdata[11];
k1[4] = steps*interdata[4];
k1[5] = steps*interdata[5];
k1[6] = steps*interdata[6];
k1[7] = steps*interdata[7];

interdata[0] = interdata[0]+0.5*k1[4];
interdata[1] = interdata[1]+0.5*k1[5];
interdata[2] = interdata[2]+0.5*k1[6];
interdata[3] = interdata[3]+0.5*k1[7];
interdata[4] = interdata[4]+0.5*k1[0];
interdata[5] = interdata[5]+0.5*k1[1];
interdata[6] = interdata[6]+0.5*k1[2];
interdata[7] = interdata[7]+0.5*k1[3];
interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);

//k2

k2[0] = steps*interdata[8];
k2[1] = steps*interdata[9];
k2[2] = steps*interdata[10];
k2[3] = steps*interdata[11];
k2[4] = steps*interdata[4];
k2[5] = steps*interdata[5];
k2[6] = steps*interdata[6];
k2[7] = steps*interdata[7];

interdata[0] = interdata[0]+0.5*k2[4];
interdata[1] = interdata[1]+0.5*k2[5];
interdata[2] = interdata[2]+0.5*k2[6];
interdata[3] = interdata[3]+0.5*k2[7];
interdata[4] = interdata[4]+0.5*k2[0];
interdata[5] = interdata[5]+0.5*k2[1];
interdata[6] = interdata[6]+0.5*k2[2];
interdata[7] = interdata[7]+0.5*k2[3];
interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);

//k3

k3[0] = steps*interdata[8];
k3[1] = steps*interdata[9];
k3[2] = steps*interdata[10];
k3[3] = steps*interdata[11];
k3[4] = steps*interdata[4];
k3[5] = steps*interdata[5];
k3[6] = steps*interdata[6];
k3[7] = steps*interdata[7];

interdata[0] = interdata[0]+k3[4];
interdata[1] = interdata[1]+k3[5];
interdata[2] = interdata[2]+k3[6];
interdata[3] = interdata[3]+k3[7];
interdata[4] = interdata[4]+k3[0];
interdata[5] = interdata[5]+k3[1];
interdata[6] = interdata[6]+k3[2];
interdata[7] = interdata[7]+k3[3];
interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);

//k4

k4[0] = steps*interdata[8];
k4[1] = steps*interdata[9];
k4[2] = steps*interdata[10];
k4[3] = steps*interdata[11];
k4[4] = steps*interdata[4];
k4[5] = steps*interdata[5];
k4[6] = steps*interdata[6];
k4[7] = steps*interdata[7];

data[0][i+1] = data[0][i]+(1.0/6.0)*(k1[0]+2*k2[0]+2*k3[0]+k4[0]);
data[1][i+1] = data[1][i]+(1.0/6.0)*(k1[1]+2*k2[1]+2*k3[1]+k4[1]);
data[2][i+1] = data[2][i]+(1.0/6.0)*(k1[2]+2*k2[2]+2*k3[2]+k4[2]);
data[3][i+1] = data[3][i]+(1.0/6.0)*(k1[3]+2*k2[3]+2*k3[3]+k4[3]);
data[4][i+1] = data[4][i]+(1.0/6.0)*(k1[4]+2*k2[4]+2*k3[4]+k4[4]);
data[5][i+1] = data[5][i]+(1.0/6.0)*(k1[5]+2*k2[5]+2*k3[5]+k4[5]);
data[6][i+1] = data[6][i]+(1.0/6.0)*(k1[6]+2*k2[6]+2*k3[6]+k4[6]);
data[7][i+1] = data[7][i]+(1.0/6.0)*(k1[7]+2*k2[7]+2*k3[7]+k4[7]);

}
//output new values
for (int i=0; i<datalength; i++){
outfile1<<i*steps<<'\t'<<data[0][i]<<'\t'<<data[1][i]<<'\t'<<data[4][i]<<'\t'<<data[5][i]<<'\t'<<data[2][i]<<'\t'<<data[3][i]<<'\t'<<data[6][i]<<'\t'<<data[7][i]<<endl;
}
outfile1.close();
}
 PhysOrg.com science news on PhysOrg.com >> City-life changes blackbird personalities, study shows>> Origins of 'The Hoff' crab revealed (w/ Video)>> Older males make better fathers: Mature male beetles work harder, care less about female infidelity
 Recognitions: Homework Help I'm no C++ expert, but I think there's something funny when void RK4_func has a variable list which consists of 'double,double,double, double, etc.) Shouldn't there be some actual parameter names in there somewhere?

 Quote by SteamKing I'm no C++ expert, but I think there's something funny when void RK4_func has a variable list which consists of 'double,double,double, double, etc.) Shouldn't there be some actual parameter names in there somewhere?
Forward declarations don't need them. The actual function is much farther down.

## Two Body problem with RK4 in C++

OP, aren't you getting errors in your compiler? What do they say? You should be, because I'm getting a ton.

Mentor
 Quote by atsquishy But my program is not compiling. Any suggestions/problems that you can see?
Don't you think it would help if you tell us the error messages that the compiler is giving you?
 You have lots of typos and some more principal problems. Always start from the top of your compiler error output. Correct simple typos (commas, semicolons, etc.) Exponentiation is different in c++, you need to import pow() function or overload ^ as such beforehand.
 Sorry sorry, yes I've made some fundamental problems. Not very used to C++, I'm rewriting it now and I'll post up the new code. I'll also be deleting the forward declaration, and making it easier to read.

Recognitions:
Science Advisor
 Quote by SteamKing I'm no C++ expert, but I think there's something funny when void RK4_func has a variable list which consists of 'double,double,double, double, etc.) Shouldn't there be some actual parameter names in there somewhere?
People often put variable names in function declarations, but strictly speaking it's better not to.

Think what could happen if a "variable name" was redefined in a header file that you included, for example. At best you would get a very confusing error message. At worst, the code would be valid, but not what you thought you had written.
 I redid it, but I still get this error. TwoBodies.cpp: In function ‘void RK4_func(double, double, double, double, double, double, double, double, double, double, double, double)’: TwoBodies.cpp:87: error: size of array ‘data’ has non-integral type ‘double’ line 87 is double data [8][datalength]; The code is Code: #include #include #include #include using namespace std; void RK4_func (double, double, double, double, double, double, double, double, double, double, double, double); int main(){ //initialise variables double time_end,time_step,mass1,x1,y1,vx1,vy1,mass2,x2,y2,vx2,vy2; //Get the value into variables cout << "Enter Simulation Duration in seconds" << endl; cin >> time_end; cout << "\n"; cout << "Enter Precision (Time step. Recommended 100 to 300 seconds )" << endl; cin >> time_step; cout << "\n"; cout << "Enter mass of body 1 with 10^20kg scale" << endl; cin >> mass1; cout << "\n"; cout << "Enter initial body 1 x position" << endl; cin >> x1; cout << "\n"; cout << "Enter initial body 1 y position" << endl; cin >> y1; cout << "\n"; cout << "Enter initial body 1 x velocity" << endl; cin >> vx1; cout << "\n"; cout << "Enter initial body 1 x velocity" << endl; cin >> vy1; cout << "\n"; cout << "Enter mass of body 2 with 10^20kg scale" << endl; cin >> mass2; cout << "\n"; cout << "Enter initial body 2 x position" << endl; cin >> x2; cout << "\n"; cout << "Enter initial body 2 y position" << endl; cin >> y2; cout << "\n"; cout << "Enter initial body 2 x velocity" << endl; cin >> vx2; cout << "\n"; cout << "Enter initial body 2 y velocity" << endl; cin >> vy2; cout << "\n"; //Pass values into RK4_func and Output data via Excel File RK4_func(time_end,time_step,mass1,x1,y1,vx1,vy1,mass2,x2,y2,vx2,vy2); return 0; } //Assume function enters data in SI units void RK4_func (double timeend, double steps, double m1, double x1, double y1, double vx1, double vy1, double m2, double x2, double y2, double vx2,double vy2) { //Excel file creation ofstream outfile1 ("RK4.xls"); outfile1<<"time"<<'\t'<<"x1"<<'\t'<<"y1"<<'\t'<<"vx1"<<'\t'<<"vy1"<<'\t'<<"x2"<<'\t'<<"y2"<<'\t'<<"vx2"<<'\t'<<"vy2"< interdata(12), k1 (8), k2 (8), k3 (8), k4 (8); // let 0~3 indicate position x1y1x2y2 4~7 indicate vx1vy1vx2vy2 8~11 indicate ax1ay1ax2ay2 //acceleration x1 y1 x2 y2 interdata[0]=data[0][i]; interdata[1]=data[1][i]; interdata[2]=data[2][i]; interdata[3]=data[3][i]; interdata[4]=data[4][i]; interdata[5]=data[5][i]; interdata[6]=data[6][i]; interdata[7]=data[7][i]; interdata[8] = G*m2/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[2][i]-data[0][i]); interdata[9] = G*m2/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[3][i]-data[1][i]); interdata[10] = G*m1/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[0][i]-data[2][i]); interdata[11] = G*m1/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[1][i]-data[3][i]); //k1 is Euler Step k1[0] = steps*interdata[8]; k1[1] = steps*interdata[9]; k1[2] = steps*interdata[10]; k1[3] = steps*interdata[11]; k1[4] = steps*interdata[4]; k1[5] = steps*interdata[5]; k1[6] = steps*interdata[6]; k1[7] = steps*interdata[7]; interdata[0] = interdata[0]+0.5*k1[4]; interdata[1] = interdata[1]+0.5*k1[5]; interdata[2] = interdata[2]+0.5*k1[6]; interdata[3] = interdata[3]+0.5*k1[7]; interdata[4] = interdata[4]+0.5*k1[0]; interdata[5] = interdata[5]+0.5*k1[1]; interdata[6] = interdata[6]+0.5*k1[2]; interdata[7] = interdata[7]+0.5*k1[3]; interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]); interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]); interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]); interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]); //k2 k2[0] = steps*interdata[8]; k2[1] = steps*interdata[9]; k2[2] = steps*interdata[10]; k2[3] = steps*interdata[11]; k2[4] = steps*interdata[4]; k2[5] = steps*interdata[5]; k2[6] = steps*interdata[6]; k2[7] = steps*interdata[7]; interdata[0] = interdata[0]+0.5*k2[4]; interdata[1] = interdata[1]+0.5*k2[5]; interdata[2] = interdata[2]+0.5*k2[6]; interdata[3] = interdata[3]+0.5*k2[7]; interdata[4] = interdata[4]+0.5*k2[0]; interdata[5] = interdata[5]+0.5*k2[1]; interdata[6] = interdata[6]+0.5*k2[2]; interdata[7] = interdata[7]+0.5*k2[3]; interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]); interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]); interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]); interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]); //k3 k3[0] = steps*interdata[8]; k3[1] = steps*interdata[9]; k3[2] = steps*interdata[10]; k3[3] = steps*interdata[11]; k3[4] = steps*interdata[4]; k3[5] = steps*interdata[5]; k3[6] = steps*interdata[6]; k3[7] = steps*interdata[7]; interdata[0] = interdata[0]+k3[4]; interdata[1] = interdata[1]+k3[5]; interdata[2] = interdata[2]+k3[6]; interdata[3] = interdata[3]+k3[7]; interdata[4] = interdata[4]+k3[0]; interdata[5] = interdata[5]+k3[1]; interdata[6] = interdata[6]+k3[2]; interdata[7] = interdata[7]+k3[3]; interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]); interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]); interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]); interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]); //k4 k4[0] = steps*interdata[8]; k4[1] = steps*interdata[9]; k4[2] = steps*interdata[10]; k4[3] = steps*interdata[11]; k4[4] = steps*interdata[4]; k4[5] = steps*interdata[5]; k4[6] = steps*interdata[6]; k4[7] = steps*interdata[7]; //Form the next position and velocity data[0][i+1] = data[0][i]+(1.0/6.0)*(k1[0]+2*k2[0]+2*k3[0]+k4[0]); data[1][i+1] = data[1][i]+(1.0/6.0)*(k1[1]+2*k2[1]+2*k3[1]+k4[1]); data[2][i+1] = data[2][i]+(1.0/6.0)*(k1[2]+2*k2[2]+2*k3[2]+k4[2]); data[3][i+1] = data[3][i]+(1.0/6.0)*(k1[3]+2*k2[3]+2*k3[3]+k4[3]); data[4][i+1] = data[4][i]+(1.0/6.0)*(k1[4]+2*k2[4]+2*k3[4]+k4[4]); data[5][i+1] = data[5][i]+(1.0/6.0)*(k1[5]+2*k2[5]+2*k3[5]+k4[5]); data[6][i+1] = data[6][i]+(1.0/6.0)*(k1[6]+2*k2[6]+2*k3[6]+k4[6]); data[7][i+1] = data[7][i]+(1.0/6.0)*(k1[7]+2*k2[7]+2*k3[7]+k4[7]); } //output new values for (int i=0; i
 You must used a fixed point type as an index into an array, you can't use a floating point value. Likewise you must have a fixed point count of objects (same thing). You can't have 3.2 elements in the array...
 so how do i specify that? do i use long int or int? EDIT: Yes. I switched the data into unsigned long int, then it compiled! but now the RK4 seems off since it doesn't give me a result.
 You have no experience in programming? It sounds like you don't even know where to begin when it comes to debugging.
 Sorry, I don't really have that much experience. It's my first course, and it's only been a few months since I started learning C++/programming. But it's compiling now. I just think the RK4 is off now, because it gives me the excel file, but the numbers are ridiculous. Is there something that I did wrong with my RK4? I mean, i'm sure there must be, but I don't see it. The RK4 function is here: Code: void RK4_func (double timeend, double steps, double m1, double x1, double y1, double vx1, double vy1, double m2, double x2, double y2, double vx2,double vy2) { //Excel file creation ofstream outfile1 ("RK4.xls"); outfile1<<"time"<<'\t'<<"x1"<<'\t'<<"y1"<<'\t'<<"vx1"<<'\t'<<"vy1"<<'\t'<<"x2"<<'\t'<<"y2"<<'\t'<<"vx2"<<'\t'<<"vy2"< interdata(12), k1 (8), k2 (8), k3 (8), k4 (8); // let 0~3 indicate position x1y1x2y2 4~7 indicate vx1vy1vx2vy2 8~11 indicate ax1ay1ax2ay2 //acceleration x1 y1 x2 y2 interdata[0]=data[0][i]; interdata[1]=data[1][i]; interdata[2]=data[2][i]; interdata[3]=data[3][i]; interdata[4]=data[4][i]; interdata[5]=data[5][i]; interdata[6]=data[6][i]; interdata[7]=data[7][i]; interdata[8] = G*m2/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[2][i]-data[0][i]); interdata[9] = G*m2/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[3][i]-data[1][i]); interdata[10] = G*m1/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[0][i]-data[2][i]); interdata[11] = G*m1/pow(pow((pow((data[0][i]-data[2][i]),2)+pow((data[1][i]-data[3][i]),2)),(1/2)),3)*(data[1][i]-data[3][i]); //k1 is Euler Step k1[0] = steps*interdata[8]; k1[1] = steps*interdata[9]; k1[2] = steps*interdata[10]; k1[3] = steps*interdata[11]; k1[4] = steps*interdata[4]; k1[5] = steps*interdata[5]; k1[6] = steps*interdata[6]; k1[7] = steps*interdata[7]; interdata[0] = interdata[0]+0.5*k1[4]; interdata[1] = interdata[1]+0.5*k1[5]; interdata[2] = interdata[2]+0.5*k1[6]; interdata[3] = interdata[3]+0.5*k1[7]; interdata[4] = interdata[4]+0.5*k1[0]; interdata[5] = interdata[5]+0.5*k1[1]; interdata[6] = interdata[6]+0.5*k1[2]; interdata[7] = interdata[7]+0.5*k1[3]; interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]); interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]); interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]); interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]); //k2 k2[0] = steps*interdata[8]; k2[1] = steps*interdata[9]; k2[2] = steps*interdata[10]; k2[3] = steps*interdata[11]; k2[4] = steps*interdata[4]; k2[5] = steps*interdata[5]; k2[6] = steps*interdata[6]; k2[7] = steps*interdata[7]; interdata[0] = interdata[0]+0.5*k2[4]; interdata[1] = interdata[1]+0.5*k2[5]; interdata[2] = interdata[2]+0.5*k2[6]; interdata[3] = interdata[3]+0.5*k2[7]; interdata[4] = interdata[4]+0.5*k2[0]; interdata[5] = interdata[5]+0.5*k2[1]; interdata[6] = interdata[6]+0.5*k2[2]; interdata[7] = interdata[7]+0.5*k2[3]; interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]); interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]); interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]); interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]); //k3 k3[0] = steps*interdata[8]; k3[1] = steps*interdata[9]; k3[2] = steps*interdata[10]; k3[3] = steps*interdata[11]; k3[4] = steps*interdata[4]; k3[5] = steps*interdata[5]; k3[6] = steps*interdata[6]; k3[7] = steps*interdata[7]; interdata[0] = interdata[0]+k3[4]; interdata[1] = interdata[1]+k3[5]; interdata[2] = interdata[2]+k3[6]; interdata[3] = interdata[3]+k3[7]; interdata[4] = interdata[4]+k3[0]; interdata[5] = interdata[5]+k3[1]; interdata[6] = interdata[6]+k3[2]; interdata[7] = interdata[7]+k3[3]; interdata[8] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[2]-interdata[0]); interdata[9] = G*m2/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[3]-interdata[1]); interdata[10] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[0]-interdata[2]); interdata[11] = G*m1/pow(pow((pow((interdata[0]-interdata[2]),2)+pow((interdata[1]-interdata[3]),2)),(1/2)),3)*(interdata[1]-interdata[3]); //k4 k4[0] = steps*interdata[8]; k4[1] = steps*interdata[9]; k4[2] = steps*interdata[10]; k4[3] = steps*interdata[11]; k4[4] = steps*interdata[4]; k4[5] = steps*interdata[5]; k4[6] = steps*interdata[6]; k4[7] = steps*interdata[7]; //Form the next position and velocity data[0][i+1] = data[0][i]+(1.0/6.0)*(k1[0]+2*k2[0]+2*k3[0]+k4[0]); data[1][i+1] = data[1][i]+(1.0/6.0)*(k1[1]+2*k2[1]+2*k3[1]+k4[1]); data[2][i+1] = data[2][i]+(1.0/6.0)*(k1[2]+2*k2[2]+2*k3[2]+k4[2]); data[3][i+1] = data[3][i]+(1.0/6.0)*(k1[3]+2*k2[3]+2*k3[3]+k4[3]); data[4][i+1] = data[4][i]+(1.0/6.0)*(k1[4]+2*k2[4]+2*k3[4]+k4[4]); data[5][i+1] = data[5][i]+(1.0/6.0)*(k1[5]+2*k2[5]+2*k3[5]+k4[5]); data[6][i+1] = data[6][i]+(1.0/6.0)*(k1[6]+2*k2[6]+2*k3[6]+k4[6]); data[7][i+1] = data[7][i]+(1.0/6.0)*(k1[7]+2*k2[7]+2*k3[7]+k4[7]); } //output new values for (int i=0; i
 Being confused is normal. Knowing how to dissect the problems are the key. To be honest, I don't understand your code. k1[0] to k1[3] should give the derivative of coordinates (velocity) times the stepsize, but you have second derivative instead (the acceleration). Same goes for k1[4] to k[7]. Shouldn't they have stepsize times the acceleration instead? Of course then you seem to flip it back to correct places in interdata, but the final calculation uses k vectors where the problem remains. Maybe it's possible to correct the problem by only changing the coefficient vector numbering in data[0][i+1] = data[0][i] + ... ... ... etc If all else fails try doing only Euler step method at first ( comment the rest out). Then there's less to think about.

 Tags c++, kepler's problem, physics, programming, rk4

 Similar discussions for: Two Body problem with RK4 in C++ Thread Forum Replies Classical Physics 2 Introductory Physics Homework 2 Special & General Relativity 2 Astrophysics 1 General Physics 9