Two Body problem with RK4 in C++

atsquishy · Dec 14, 2012

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();
}

SteamKing · Dec 14, 2012

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?

justsomeguy · Dec 14, 2012

SteamKing said:

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.

justsomeguy · Dec 14, 2012

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

jtbell · Dec 14, 2012

atsquishy said:

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?

Martin_L · Dec 15, 2012

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.

atsquishy · Dec 15, 2012

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.

AlephZero · Dec 15, 2012

SteamKing said:

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.

atsquishy · Dec 15, 2012

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

justsomeguy · Dec 15, 2012

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...

atsquishy · Dec 15, 2012

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.

justsomeguy · Dec 15, 2012

You have no experience in programming? It sounds like you don't even know where to begin when it comes to debugging.

atsquishy · Dec 15, 2012

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();
}

Martin_L · Dec 19, 2012

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] + ...
...
...
etc

If all else fails try doing only Euler step method at first ( comment the rest out). Then there's less to think about.

Two Body problem with RK4 in C++

1. What is the Two Body problem in physics?

2. What is RK4 in C++?

3. How does RK4 work in solving the Two Body problem?

4. What is the advantage of using RK4 in C++ for the Two Body problem?

5. How can I implement RK4 in C++ for the Two Body problem?

Similar threads

Hot Threads

Recent Insights