Two Body problem with RK4 in C++


by atsquishy
Tags: c++, kepler's problem, physics, programming, rk4
atsquishy
atsquishy is offline
#1
Dec14-12, 06:09 PM
P: 5
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?

#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();
}
Phys.Org News Partner Science news on Phys.org
Simplicity is key to co-operative robots
Chemical vapor deposition used to grow atomic layer materials on top of each other
Earliest ancestor of land herbivores discovered
SteamKing
SteamKing is offline
#2
Dec14-12, 07:09 PM
HW Helper
Thanks
P: 5,529
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
justsomeguy is offline
#3
Dec14-12, 07:14 PM
P: 166
Quote Quote by SteamKing View Post
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
justsomeguy is offline
#4
Dec14-12, 07:15 PM
P: 166

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.
jtbell
jtbell is offline
#5
Dec14-12, 08:12 PM
Mentor
jtbell's Avatar
P: 11,221
Quote Quote by atsquishy View Post
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
Martin_L is offline
#6
Dec15-12, 03:35 AM
P: 5
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
atsquishy is offline
#7
Dec15-12, 05:48 AM
P: 5
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
AlephZero is offline
#8
Dec15-12, 11:48 AM
Engineering
Sci Advisor
HW Helper
Thanks
P: 6,341
Quote Quote by SteamKing View Post
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
atsquishy is offline
#9
Dec15-12, 05:07 PM
P: 5
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
#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
justsomeguy is offline
#10
Dec15-12, 05:21 PM
P: 166
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
atsquishy is offline
#11
Dec15-12, 07:25 PM
P: 5
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
justsomeguy is offline
#12
Dec15-12, 07:45 PM
P: 166
You have no experience in programming? It sounds like you don't even know where to begin when it comes to debugging.
atsquishy
atsquishy is offline
#13
Dec15-12, 08:15 PM
P: 5
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:
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
Martin_L is offline
#14
Dec19-12, 05:04 AM
P: 5
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.


Register to reply

Related Discussions
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