| New Reply |
Two Body problem with RK4 in C++ |
Share Thread | Thread Tools |
| Dec14-12, 06:09 PM | #1 |
|
|
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();
}
|
| Dec14-12, 07:09 PM | #2 |
|
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? |
| Dec14-12, 07:14 PM | #3 |
|
|
|
|
| Dec14-12, 07:15 PM | #4 |
|
|
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.
|
| Dec14-12, 08:12 PM | #5 |
|
Mentor
|
|
|
| Dec15-12, 03:35 AM | #6 |
|
|
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. |
| Dec15-12, 05:48 AM | #7 |
|
|
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.
|
| Dec15-12, 11:48 AM | #8 |
Recognitions:
Science Advisor |
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. |
|
| Dec15-12, 05:07 PM | #9 |
|
|
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 <iostream>
#include <vector>
#include <fstream>
#include <math.h>
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"<<endl;
outfile1<<0<<'\t'<<x1<<'\t'<<y1<<'\t'<<vx1<<'\t'<<vy1<<'\t'<<x2<<'\t'<<y2<<'\t'<<vx2<<'\t'<<vy2<<endl;
//Declare constants and array
const double G = 6.67*pow(10,-11);
double datalength = 1+(timeend/steps);
double data[8][datalength];
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/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<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();
}
// Display via GnuPlot
|
| Dec15-12, 05:21 PM | #10 |
|
|
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...
|
| Dec15-12, 07:25 PM | #11 |
|
|
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. |
| Dec15-12, 07:45 PM | #12 |
|
|
You have no experience in programming? It sounds like you don't even know where to begin when it comes to debugging.
|
| Dec15-12, 08:15 PM | #13 |
|
|
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"<<endl;
outfile1<<0<<'\t'<<x1<<'\t'<<y1<<'\t'<<vx1<<'\t'<<vy1<<'\t'<<x2<<'\t'<<y2<<'\t'<<vx2<<'\t'<<vy2<<endl;
//Declare constants and array
const double G = 6.67*pow(10,-11);
unsigned long int datalength = 1+(timeend/steps);
double data[8][datalength];
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/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<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();
}
|
| Dec19-12, 05:04 AM | #14 |
|
|
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. |
| New Reply |
| Tags |
| c++, kepler's problem, physics, programming, rk4 |
| Thread Tools | |
Similar Threads for: Two Body problem with RK4 in C++
|
||||
| Thread | Forum | Replies | ||
| Is the energy conserved FOR EACH BODY in a two-body central force problem? | Classical Physics | 2 | ||
| Not sure how to do 2 body problem | Introductory Physics Homework | 2 | ||
| two body problem in GR | Special & General Relativity | 2 | ||
| Three body problem | Astrophysics | 1 | ||
| N-body problem | General Physics | 9 | ||
Physorg.com Science News Partner
Quote by SteamKing
