I am working on a famous classical physics problem: - the three body problem. I have made the follwoing program : -
#include <math.h>
#include <iostream>
using namespace std;
class System_variables
{
private:
// These are the system variables:- position, velocity and accleration
long double body_1_pos[3];
long double body_2_pos[3];
long double body_3_pos[3];
long double body_1_v[3];
long double body_2_v[3];
long double body_3_v[3];
long double acc_body_1[3];
long double acc_body_2[3];
long double acc_body_3[3];
// Needed constants
const long double body_mass = 6 * pow(10, 24);
double time_step;
const long double G = 6.67 * pow(10, -11);
public:
System_variables(long double body_1_pos_i[3], long double body_2_pos_i[3], long double body_3_pos_i[3], long double body_1_v_i[3], long double body_2_v_i[3], long double body_3_v_i[3])
{
// Initiates the system variables, position and velocity.
for (int m = 0; m<=2; m++)
{
body_1_pos[m] = body_1_pos_i[m];
body_2_pos[m] = body_2_pos_i[m];
body_3_pos[m] = body_3_pos_i[m];
body_1_v[m] = body_1_v_i[m];
body_2_v[m] = body_2_v_i[m];
body_3_v[m] = body_3_v_i[m];
}
}
void init_time_step(double step)
{
// Initializes time step.
time_step = step;
}
void calculate_acc()
{
// Calculates the accleration on the bodies using newtons equations.
for (int m=0; m<=2; m++)
{
acc_body_1[m] = G * body_mass/pow((body_2_pos[m] - body_1_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_1_pos[m]), 2);
acc_body_2[m] = G * body_mass/pow((body_1_pos[m] - body_2_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_2_pos[m]), 2);
acc_body_3[m] = G * body_mass/pow((body_1_pos[m] - body_3_pos[m]), 2) + G * body_mass/pow((body_2_pos[m] - body_3_pos[m]), 2);
};
return;
}
void apply_acc_to_v()
{
// Applys the accleration to the velocities of the bodies.
for (int m=0; m<=2; m++)
{
body_1_v[m] += acc_body_1[m] * time_step;
body_2_v[m] += acc_body_2[m] * time_step;
body_3_v[m] += acc_body_3[m] * time_step;
}
return;
}
void apply_v_to_pos()
{
// Applys the velocities to the positions of the bodies.
for (int m=0; m<=2; m++)
{
body_1_pos[m] += body_1_v[m] * time_step;
body_2_pos[m] += body_2_v[m] * time_step;
body_3_pos[m] += body_3_v[m] * time_step;
}
return;
}
long double get_body_1_pos(int index)
{
// Returns the position of body 1 based on index.
return body_1_pos[index];
}
long double get_body_2_pos(int index)
{
// Returns the position of body 2 based on index.
return body_2_pos[index];
}
long double get_body_3_pos(int index)
{
// Returns the position of body 3 based on index.
return body_3_pos[index];
}
long double get_body_1_v(int index)
{
// Returns the velocity of body 1 based on index.
return body_1_v[index];
}
long double get_body_2_v(int index)
{
// Returns the velocity of body 2 based on index.
return body_2_v[index];
}
long double get_body_3_v(int index)
{
// Returns the velocity of body 3 based on index.
return body_3_v[index];
}
};
int main()
{
long double array_1[3] = {100,100,100};
long double array_2[3] ={10000000000, 10000000000, 10000000000};
long double array_3[3] = {10000000000, 10000000000, -10000000000};
long double array_0[3] = {0,0,0};
System_variables system(array_1, array_2, array_3, array_0, array_0, array_0);
system.init_time_step(100);
for (int time = 0; time <= 10; time++)
{
cout << "x position of body 1: " << system.get_body_1_pos(0) << endl
<< "y position of body 1: " << system.get_body_1_pos(1) << endl
<< "z position of body 1: " << system.get_body_1_pos(2) << endl << endl;
cout << "x position of body 2: " << system.get_body_2_pos(0) << endl
<< "y position of body 2: " << system.get_body_2_pos(1) << endl
<< "z position of body 2: " << system.get_body_2_pos(2) << endl << endl;
cout << "x position of body 3: " << system.get_body_3_pos(0) << endl
<< "y position of body 3: " << system.get_body_3_pos(1) << endl
<< "z position of body 3: " << system.get_body_3_pos(2) << endl << endl;
system.calculate_acc();
system.apply_acc_to_v();
system.apply_v_to_pos();
}
}
A point to note is that instead of using the classic double data type, I have used long double, because when I used the double data type, after a while inf was being outputted, so I used the long double thinking that it would store decimals to large precision. But, it didn't work. I tried also float (though it has less precision), but the output is always the same : -
x position of body 1: 100
y position of body 1: 100
z position of body 1: 100
x position of body 2: 1e+10
y position of body 2: 1e+10
z position of body 2: 1e+10
x position of body 3: 1e+10
y position of body 3: 1e+10
z position of body 3: -1e+10
x position of body 1: 100.08
y position of body 1: 100.08
z position of body 1: 100.08
x position of body 2: inf
y position of body 2: inf
z position of body 2: 1e+10
x position of body 3: inf
y position of body 3: inf
z position of body 3: -1e+10
x position of body 1: 100.16
y position of body 1: 100.16
z position of body 1: 100.24
x position of body 2: -nan(ind)
y position of body 2: -nan(ind)
z position of body 2: 1e+10
x position of body 3: -nan(ind)
y position of body 3: -nan(ind)
z position of body 3: -1e+10
x position of body 1: -nan(ind)
y position of body 1: -nan(ind)
z position of body 1: 100.48
x position of body 2: -nan(ind)
y position of body 2: -nan(ind)
z position of body 2: 1e+10
x position of body 3: -nan(ind)
y position of body 3: -nan(ind)
z position of body 3: -1e+10
I am thinking of using pointers, but I doubt if it is a right idea. What do you think?
Your calculations are not numerically stable.
The exact place where things start to break first in your program is the acceleration calculation. For instance:
acc_body_2[m] = G * body_mass/pow((body_1_pos[m] - body_2_pos[m]), 2) +
G * body_mass/pow((body_3_pos[m] - body_2_pos[m]), 2);
Body_3 and Body_2 contain very big numbers of similar magnitude. Their difference will be a very small number. Squaring that small number makes it even smaller. You then divide the mass by that very small number, creating a very,very big number. This is not good.
In your case the numbers for body_3 and body_2 are even identical, so you're really dividing by 0 here, which sends you straight to inf. But even with numbers that are only somewhat equal, you will likely run into problems here.
Note that this is not primarily a problem of the limited precision of the used floating point data type (although the problem can get mitigated to some extent by adding more precision), but primarily a problem of the way you are structuring out computation. There is an entire discipline of computer science called numerical computing dedicated to deal with the subtle issues that can arise with computations like this.
I'd suggest you start by doing some research into numerical integration, which is really what you are trying to do here, and look at some of the common pitfalls like the catastrophic cancellation problem.