Two Body problem with RK4 in C++

  • C/++/#
  • Thread starter atsquishy
  • Start date
  • #1
5
0

Main Question or Discussion Point

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

Answers and Replies

  • #2
SteamKing
Staff Emeritus
Science Advisor
Homework Helper
12,796
1,667
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?
 
  • #3
153
0
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.
 
  • #4
153
0
OP, aren't you getting errors in your compiler? What do they say? You should be, because I'm getting a ton.
 
  • #5
jtbell
Mentor
15,581
3,560
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? :confused:
 
  • #6
6
0
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.
 
  • #7
5
0
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.
 
  • #8
AlephZero
Science Advisor
Homework Helper
6,994
291
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.
 
  • #9
5
0
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
 
Last edited:
  • #10
153
0
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...
 
  • #11
5
0
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.
 
  • #12
153
0
You have no experience in programming? It sounds like you don't even know where to begin when it comes to debugging.
 
  • #13
5
0
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();
}
 
Last edited:
  • #14
6
0
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.
 

Related Threads on Two Body problem with RK4 in C++

Replies
11
Views
7K
Replies
3
Views
6K
  • Last Post
Replies
1
Views
3K
  • Last Post
Replies
1
Views
5K
Replies
7
Views
948
  • Last Post
Replies
6
Views
1K
  • Last Post
Replies
1
Views
1K
Replies
3
Views
2K
Replies
2
Views
1K
  • Last Post
Replies
18
Views
6K
Top