EnergumeniEnergy of the system with the proposed solution:
[
]
Energumeni of the system with the proposed solution:
[]
Energy of the system with the proposed solution:
[]
Edit 2;
Energumeni of the system with the proposed solution:
[]
Edit 2;
Energumeni of the system with the proposed solution:
[]
#include <iostream>
#include <cmath>
#include <array>
#include <fstream>
#include "integrators.h"
#include <string>
#include <map>
// Value-Defintions of the different String values
enum StringValue { evNotDefined,
evStringValue1,
evStringValue2
};
// Map to associate the strings with the enum values
static std::map<std::string, StringValue> s_mapStringValues;
void Initialize(){
s_mapStringValues["euler"] = evStringValue1;
s_mapStringValues["rk4"] = evStringValue2;
}
static constexpr int DIM = 4;
static constexpr double G = 10;
static constexpr int N_BODIES = 3;
static constexpr int N_STEPS = 300000;
double distance(std::array<double, 3> r1, std::array<double, 3> r2){
return sqrt(pow(r1[0]-r2[0],2)+pow(r1[1]-r2[1],2)+pow(r1[2]-r2[2],2));
}
class Planet{
public:
double m;
std::array <double, 3> x;
std::array <double, 3> v;
std::array <double, 3> a;
double energy;
Planet (double mass, double x_position, double y_position, double z_position, double x_velocity, double y_velocity, double z_velocity) {
m = mass;
x[0] = x_position;
x[1] = y_position;
x[2] = z_position;
v[0] = x_velocity;
v[1] = y_velocity;
v[2] = z_velocity;
};
void computeKineticEnergy(){
energy = 0.5 * m * (pow(v[0],2)+pow(v[1],2)+pow(v[2],2));
}
};
double computePotentialEnergy(Planet planet1, Planet planet2, Planet planet3){
return -1 * G * ( planet3.m * planet1.m / distance(planet3.x, planet1.x) + planet3.m * planet2.m / distance(planet3.x, planet2.x) + planet1.m * planet2.m / distance(planet1.x, planet2.x));
}
double acceleration(Planet A, Planet B, Planet C, int axe){
// Compute the acceleration of the body C, specifying the axis.
double mass_A = A.m;
double mass_B = B.m;
double posx_A = A.x[0];
double posx_B = B.x[0];
double posx_C = C.x[0];
double posy_A = A.x[1];
double posy_B = B.x[1];
double posy_C = C.x[1];
double posz_A = A.x[2];
double posz_B = B.x[2];
double posz_C = C.x[2];
if (axe == 0){
return (-1 * G * (mass_A * (posx_C-posx_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posx_C-posx_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (axe == 1){
return (-1 * G * (mass_A * (posy_C-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posy_C-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (axe == 2){
return (-1 * G * (mass_A * (posz_C-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posz_C-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}
}
double F(double x, double v, double t, Planet A, Planet B, Planet C, int j ){
// Function to integrate via Runge-Kutta.
double mass_A = A.m;
double mass_B = B.m;
double posx_A = A.x[0];
double posx_B = B.x[0];
double posx_C = C.x[0];
double posy_A = A.x[1];
double posy_B = B.x[1];
double posy_C = C.x[1];
double posz_A = A.x[2];
double posz_B = B.x[2];
double posz_C = C.x[2];
if (j == 0){
return (-1 * G * (mass_A * (x-posx_A) / pow(sqrt(pow(x-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posx_B) / pow(sqrt(pow(x-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 1) {
return (-1 * G * (mass_A * (x-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(x-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(x-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 2) {
return (-1 * G * (mass_A * (x-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(x-posz_A,2)), 3) + mass_B * (x-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(x-posz_B,2)), 3)));
}
}
int main(int argc, char** argv){
double h = 0.0005;
Planet A(100, -10, 10, -11, -3, 0, 0);
Planet B(100, 0, 0, 0, 3, 0, 0);
Planet C(0, 10, 14, 12, 3, 0, 0);
double x_A[DIM][3];
double x_B[DIM][3];
double x_C[DIM][3];
double vx_A;
double vy_A;
double vz_A;
double vx_B;
double vy_B;
double vz_B;
double vx_C;
double vy_C;
double vz_C;
double mass_A = A.m;
double mass_B = B.m;
double mass_C = C.m;
x_A[0][0] = A.x[0];
x_B[0][0] = B.x[0];
x_C[0][0] = C.x[0];
vx_A = A.v[0];
vx_B = B.v[0];
vx_C = C.v[0];
x_A[1][0] = A.x[1];
x_B[1][0] = B.x[1];
x_C[1][0] = C.x[1];
vy_A = A.v[1];
vy_B = B.v[1];
vy_C = C.v[1];
x_A[2][0] = A.x[2];
x_B[2][0] = B.x[2];
x_C[2][0] = C.x[2];
vz_A = A.v[2];
vz_B = B.v[2];
vz_C = C.v[2];
Initialize();
std::ofstream file_energy("Total_energy_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_A("positions_A_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_B("positions_B_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_C("positions_C_" + std::string(argv[1]) + ".csv");
output_file_A<<"x;y;z"<<std::endl;
output_file_B<<"x;y;z"<<std::endl;
output_file_C<<"x;y;z"<<std::endl;
file_energy<<"k;p"<<std::endl;
double m1;
double k1;
double m2;
double k2;
double m3;
double k3;
double m4;
double k4;
std::array<double,3> vA = A.v; //condizione iniziale velocita
std::array<double,3> xA = A.x;
std::array<double,3> vB = B.v; //condizione iniziale velocita
std::array<double,3> xB = B.x;
std::array<double,3> vC = C.v; //condizione iniziale velocita
std::array<double,3> xC = C.x;
double t;
std::array<double,3> cm = computeCM(A, B, C);
if (argc>=2){
switch (s_mapStringValues[argv[1]]){
case evStringValue1:
// ==========================================================
// EULER
// ==========================================================
for (int i=0; i<N_STEPS-1; i++){
output_file_A << A.x[0] << ";" << A.x[1] << ";" << A.x[2]<< std::endl;
output_file_B << B.x[0] << ";" << B.x[1] << ";" << B.x[2]<< std::endl;
output_file_C << C.x[0] << ";" << C.x[1] << ";" << C.x[2]<< std::endl;
for(int j=0; j<DIM-1; j++){
A.a[j] = acceleration(B, C, A, j);
B.a[j] = acceleration(A, C, B, j);
C.a[j] = acceleration(B, A, C, j);
A.v[j] += A.a[j] * h;
B.v[j] += B.a[j] * h;
C.v[j] += C.a[j] * h;
x_A[j][0] = x_A[j][0] + A.v[j] * h;
x_B[j][0] = x_B[j][0] + B.v[j] * h;
x_C[j][0] = x_C[j][0] + C.v[j] * h;
}
for (int j=0;j<DIM-1;j++){
A.x[j] = x_A[j][0];
B.x[j] = x_B[j][0];
C.x[j] = x_C[j][0];
}
A.computeKineticEnergy();
B.computeKineticEnergy();
C.computeKineticEnergy();
cm = computeCM(A,B,C);
file_energy<<A.energy + B.energy + C.energy<<";"<< computePotentialEnergy(A, B, C)<<std::endl;
}
break;
case evStringValue2:{
// ==========================================================
// RUNGE KUTTA 4
// ==========================================================
for(int i=0; i<N_STEPS-1; i++){
for(int j=0; j<DIM-1; j++){
// body A
m1 = h*vA[j];
k1 = h*F(xA[j], vA[j], t, C, B, A, j);
m2 = h*(vA[j] + 0.5*k1);
k2 = h*F(xA[j]+0.5*m1, vA[j]+0.5*k1, t+0.5*h, C, B, A, j);
m3 = h*(vA[j] + 0.5*k2);
k3 = h*F(xA[j]+0.5*m2, vA[j]+0.5*k2, t+0.5*h, C, B, A, j);
m4 = h*(vA[j] + k3);
k4 = h*F(xA[j]+m3, vA[j]+k3, t+h, C, B, A, j);
xA[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vA[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1; j++){
// body B
m1 = h*vB[j];
k1 = h*F(xB[j], vB[j], t, A, C, B, j); //(x, v, t)
m2 = h*(vB[j] + 0.5*k1);
k2 = h*F(xB[j]+0.5*m1, vB[j]+0.5*k1, t+0.5*h, A, C, B, j);
m3 = h*(vB[j] + 0.5*k2);
k3 = h*F(xB[j]+0.5*m2, vB[j]+0.5*k2, t+0.5*h, A, C, B, j);
m4 = h*(vB[j] + k3);
k4 = h*F(xB[j]+m3, vB[j]+k3, t+h, A, C, B, j);
xB[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vB[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1; j++){
// body C
m1 = h*vC[j];
k1 = h*F(xC[j], vC[j], t, A, B, C, j); //(x, v, t)
m2 = h*(vC[j] + 0.5*k1);
k2 = h*F(xC[j]+0.5*m1, vC[j]+0.5*k1, t+0.5*h, A, B, C, j);
m3 = h*(vC[j] + 0.5*k2);
k3 = h*F(xC[j]+0.5*m2, vC[j]+0.5*k2, t+0.5*h, A, B, C, j);
m4 = h*(vC[j] + k3);
k4 = h*F(xC[j]+m3, vC[j]+k3, t+h, A, B, C, j);
xC[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vC[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1;j++){
A.v[j] = vA[j];
B.v[j] = vB[j];
C.v[j] = vC[j];
A.x[j] = xA[j];
B.x[j] = xB[j];
C.x[j] = xC[j];
}
output_file_A << A.x[0] << ";" << A.x[1] << ";" << A.x[2]<< std::endl;
output_file_B << B.x[0] << ";" << B.x[1] << ";" << B.x[2]<< std::endl;
output_file_C << C.x[0] << ";" << C.x[1] << ";" << C.x[2]<< std::endl;
A.computeKineticEnergy();
B.computeKineticEnergy();
C.computeKineticEnergy();
file_energy<<A.energy + B.energy + C.energy <<";"<< computePotentialEnergy(A, B, C)<<std::endl;
}
}
break;
default:
std::cout<<"Insert an argument between: euler, or rk4"<<std::endl;
return 0;
}
}else{
std::cout<<"Insert an argument between: euler, or rk4"<<std::endl;
return 0;
}
//----------------------------------------------------------------
output_file_A.close();
output_file_B.close();
output_file_C.close();
file_energy.close();
return 0;
}
#include <iostream>
#include <cmath>
#include <array>
#include <fstream>
#include "integrators.h"
#include <string>
#include <map>
// Value-Defintions of the different String values
enum StringValue { evNotDefined,
evStringValue1,
evStringValue2
};
// Map to associate the strings with the enum values
static std::map<std::string, StringValue> s_mapStringValues;
void Initialize(){
s_mapStringValues["euler"] = evStringValue1;
s_mapStringValues["rk4"] = evStringValue2;
}
static constexpr int DIM = 4;
static constexpr double G = 10;
static constexpr int N_BODIES = 3;
static constexpr int N_STEPS = 300000;
double distance(std::array<double, 3> r1, std::array<double, 3> r2){
return sqrt(pow(r1[0]-r2[0],2)+pow(r1[1]-r2[1],2)+pow(r1[2]-r2[2],2));
}
class Planet{
public:
double m;
std::array <double, 3> x;
std::array <double, 3> v;
std::array <double, 3> a;
double energy;
Planet (double mass, double x_position, double y_position, double z_position, double x_velocity, double y_velocity, double z_velocity) {
m = mass;
x[0] = x_position;
x[1] = y_position;
x[2] = z_position;
v[0] = x_velocity;
v[1] = y_velocity;
v[2] = z_velocity;
};
void computeKineticEnergy(){
energy = 0.5 * m * (pow(v[0],2)+pow(v[1],2)+pow(v[2],2));
}
};
double computePotentialEnergy(Planet planet1, Planet planet2, Planet planet3){
return -1 * G * ( planet3.m * planet1.m / distance(planet3.x, planet1.x) + planet3.m * planet2.m / distance(planet3.x, planet2.x) + planet1.m * planet2.m / distance(planet1.x, planet2.x));
}
double acceleration(Planet A, Planet B, Planet C, int axe){
// Compute the acceleration of the body C, specifying the axis.
double mass_A = A.m;
double mass_B = B.m;
double posx_A = A.x[0];
double posx_B = B.x[0];
double posx_C = C.x[0];
double posy_A = A.x[1];
double posy_B = B.x[1];
double posy_C = C.x[1];
double posz_A = A.x[2];
double posz_B = B.x[2];
double posz_C = C.x[2];
if (axe == 0){
return (-1 * G * (mass_A * (posx_C-posx_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posx_C-posx_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (axe == 1){
return (-1 * G * (mass_A * (posy_C-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posy_C-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (axe == 2){
return (-1 * G * (mass_A * (posz_C-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posz_C-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}
}
double F(double x, double v, double t, Planet A, Planet B, Planet C, int j ){
// Function to integrate via Runge-Kutta.
double mass_A = A.m;
double mass_B = B.m;
double posx_A = A.x[0];
double posx_B = B.x[0];
double posx_C = C.x[0];
double posy_A = A.x[1];
double posy_B = B.x[1];
double posy_C = C.x[1];
double posz_A = A.x[2];
double posz_B = B.x[2];
double posz_C = C.x[2];
if (j == 0){
return (-1 * G * (mass_A * (x-posx_A) / pow(sqrt(pow(x-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posx_B) / pow(sqrt(pow(x-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 1) {
return (-1 * G * (mass_A * (x-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(x-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(x-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 2) {
return (-1 * G * (mass_A * (x-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(x-posz_A,2)), 3) + mass_B * (x-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(x-posz_B,2)), 3)));
}
}
int main(int argc, char** argv){
double h = 0.0005;
Planet A(100, -10, 10, -11, -3, 0, 0);
Planet B(100, 0, 0, 0, 3, 0, 0);
Planet C(0, 10, 14, 12, 3, 0, 0);
double x_A[DIM][3];
double x_B[DIM][3];
double x_C[DIM][3];
double vx_A;
double vy_A;
double vz_A;
double vx_B;
double vy_B;
double vz_B;
double vx_C;
double vy_C;
double vz_C;
double mass_A = A.m;
double mass_B = B.m;
double mass_C = C.m;
x_A[0][0] = A.x[0];
x_B[0][0] = B.x[0];
x_C[0][0] = C.x[0];
vx_A = A.v[0];
vx_B = B.v[0];
vx_C = C.v[0];
x_A[1][0] = A.x[1];
x_B[1][0] = B.x[1];
x_C[1][0] = C.x[1];
vy_A = A.v[1];
vy_B = B.v[1];
vy_C = C.v[1];
x_A[2][0] = A.x[2];
x_B[2][0] = B.x[2];
x_C[2][0] = C.x[2];
vz_A = A.v[2];
vz_B = B.v[2];
vz_C = C.v[2];
Initialize();
std::ofstream file_energy("Total_energy_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_A("positions_A_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_B("positions_B_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_C("positions_C_" + std::string(argv[1]) + ".csv");
output_file_A<<"x;y;z"<<std::endl;
output_file_B<<"x;y;z"<<std::endl;
output_file_C<<"x;y;z"<<std::endl;
file_energy<<"k;p"<<std::endl;
double m1;
double k1;
double m2;
double k2;
double m3;
double k3;
double m4;
double k4;
std::array<double,3> vA = A.v; //condizione iniziale velocita
std::array<double,3> xA = A.x;
std::array<double,3> vB = B.v; //condizione iniziale velocita
std::array<double,3> xB = B.x;
std::array<double,3> vC = C.v; //condizione iniziale velocita
std::array<double,3> xC = C.x;
double t;
std::array<double,3> cm = computeCM(A, B, C);
if (argc>=2){
switch (s_mapStringValues[argv[1]]){
case evStringValue1:
// ==========================================================
// EULER
// ==========================================================
for (int i=0; i<N_STEPS-1; i++){
output_file_A << A.x[0] << ";" << A.x[1] << ";" << A.x[2]<< std::endl;
output_file_B << B.x[0] << ";" << B.x[1] << ";" << B.x[2]<< std::endl;
output_file_C << C.x[0] << ";" << C.x[1] << ";" << C.x[2]<< std::endl;
for(int j=0; j<DIM-1; j++){
A.a[j] = acceleration(B, C, A, j);
B.a[j] = acceleration(A, C, B, j);
C.a[j] = acceleration(B, A, C, j);
A.v[j] += A.a[j] * h;
B.v[j] += B.a[j] * h;
C.v[j] += C.a[j] * h;
x_A[j][0] = x_A[j][0] + A.v[j] * h;
x_B[j][0] = x_B[j][0] + B.v[j] * h;
x_C[j][0] = x_C[j][0] + C.v[j] * h;
}
for (int j=0;j<DIM-1;j++){
A.x[j] = x_A[j][0];
B.x[j] = x_B[j][0];
C.x[j] = x_C[j][0];
}
A.computeKineticEnergy();
B.computeKineticEnergy();
C.computeKineticEnergy();
cm = computeCM(A,B,C);
file_energy<<A.energy + B.energy + C.energy<<";"<< computePotentialEnergy(A, B, C)<<std::endl;
}
break;
case evStringValue2:{
// ==========================================================
// RUNGE KUTTA 4
// ==========================================================
for(int i=0; i<N_STEPS-1; i++){
for(int j=0; j<DIM-1; j++){
// body A
m1 = h*vA[j];
k1 = h*F(xA[j], vA[j], t, C, B, A, j);
m2 = h*(vA[j] + 0.5*k1);
k2 = h*F(xA[j]+0.5*m1, vA[j]+0.5*k1, t+0.5*h, C, B, A, j);
m3 = h*(vA[j] + 0.5*k2);
k3 = h*F(xA[j]+0.5*m2, vA[j]+0.5*k2, t+0.5*h, C, B, A, j);
m4 = h*(vA[j] + k3);
k4 = h*F(xA[j]+m3, vA[j]+k3, t+h, C, B, A, j);
xA[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vA[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1; j++){
// body B
m1 = h*vB[j];
k1 = h*F(xB[j], vB[j], t, A, C, B, j); //(x, v, t)
m2 = h*(vB[j] + 0.5*k1);
k2 = h*F(xB[j]+0.5*m1, vB[j]+0.5*k1, t+0.5*h, A, C, B, j);
m3 = h*(vB[j] + 0.5*k2);
k3 = h*F(xB[j]+0.5*m2, vB[j]+0.5*k2, t+0.5*h, A, C, B, j);
m4 = h*(vB[j] + k3);
k4 = h*F(xB[j]+m3, vB[j]+k3, t+h, A, C, B, j);
xB[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vB[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1; j++){
// body C
m1 = h*vC[j];
k1 = h*F(xC[j], vC[j], t, A, B, C, j); //(x, v, t)
m2 = h*(vC[j] + 0.5*k1);
k2 = h*F(xC[j]+0.5*m1, vC[j]+0.5*k1, t+0.5*h, A, B, C, j);
m3 = h*(vC[j] + 0.5*k2);
k3 = h*F(xC[j]+0.5*m2, vC[j]+0.5*k2, t+0.5*h, A, B, C, j);
m4 = h*(vC[j] + k3);
k4 = h*F(xC[j]+m3, vC[j]+k3, t+h, A, B, C, j);
xC[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vC[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1;j++){
A.v[j] = vA[j];
B.v[j] = vB[j];
C.v[j] = vC[j];
A.x[j] = xA[j];
B.x[j] = xB[j];
C.x[j] = xC[j];
}
output_file_A << A.x[0] << ";" << A.x[1] << ";" << A.x[2]<< std::endl;
output_file_B << B.x[0] << ";" << B.x[1] << ";" << B.x[2]<< std::endl;
output_file_C << C.x[0] << ";" << C.x[1] << ";" << C.x[2]<< std::endl;
A.computeKineticEnergy();
B.computeKineticEnergy();
C.computeKineticEnergy();
file_energy<<A.energy + B.energy + C.energy <<";"<< computePotentialEnergy(A, B, C)<<std::endl;
}
}
break;
default:
std::cout<<"Insert an argument between: euler, or rk4"<<std::endl;
return 0;
}
}else{
std::cout<<"Insert an argument between: euler, or rk4"<<std::endl;
return 0;
}
//----------------------------------------------------------------
output_file_A.close();
output_file_B.close();
output_file_C.close();
file_energy.close();
return 0;
}
edit:
I think I have understood the problem. I have modified the code as follow, but now the energy decreases a lot and the trajectories become more and more little.
for(int i=0; i<N_STEPS-1; i++){
for(int j=0; j<DIM-1; j++){
// body A
m1[0][j] = h*vA[j];
k1[0][j] = h*F(xA[j], vA[j], t, C, B, A, j);
// body B
m1[1][j] = h*vA[j];
k1[1][j] = h*F(xA[j], vA[j], t, C, A, B, j);
// body C
m1[2][j] = h*vA[j];
k1[2][j] = h*F(xA[j], vA[j], t, A, B, C, j);
}
for(int j=0; j<DIM-1;j++){
A.v[j] += 0.5 * k1[0][j];
B.v[j] += 0.5 * k1[1][j];
C.v[j] += 0.5 * k1[2][j];
A.x[j] += 0.5 * m1[0][j];
B.x[j] += 0.5 * m1[1][j];
C.x[j] += 0.5 * m1[2][j];
}
for(int j=0; j<DIM-1; j++){
//Body A
m2[0][j] = h*(A.v[j]);
k2[0][j] = h*F(A.x[j], A.v[j], t+0.5*h, C, B, A, j);
//Body B
m2[1][j] = h*(B.v[j]);
k2[1][j] = h*F(B.x[j], B.v[j], t+0.5*h, C, A, B, j);
// Body C
m2[2][j] = h*(C.v[j]);
k2[2][j] = h*F(C.x[j], C.v[j], t+0.5*h, A, B, C, j);
}
for(int j=0; j<DIM-1;j++){
A.v[j] += 0.5 * k2[0][j];
B.v[j] += 0.5 * k2[1][j];
C.v[j] += 0.5 * k2[2][j];
A.x[j] += 0.5 * m2[0][j];
B.x[j] += 0.5 * m2[1][j];
C.x[j] += 0.5 * m2[2][j];
}
for(int j=0; j<DIM-1; j++){
//Body A
m3[0][j] = h*(A.v[j]);
k3[0][j] = h*F(A.x[j], A.v[j], t+0.5*h, C, B, A, j);
//Body B
m3[1][j] = h*(B.v[j]);
k3[1][j] = h*F(B.x[j], B.v[j], t+0.5*h, C, A, B, j);
// Body C
m3[2][j] = h*(C.v[j]);
k3[2][j] = h*F(C.x[j], C.v[j], t+0.5*h, A, B, C, j);
}
for(int j=0; j<DIM-1;j++){
A.v[j] += k3[0][j];
B.v[j] += k3[1][j];
C.v[j] += k3[2][j];
A.x[j] += m3[0][j];
B.x[j] += m3[1][j];
C.x[j] += m3[2][j];
}
for(int j=0; j<DIM-1; j++){
//Body A
m4[0][j] = h*(A.v[j]);
k4[0][j] = h*F(A.x[j], A.v[j], t + h, C, B, A, j);
//Body B
m4[1][j] = h*(B.v[j]);
k4[1][j] = h*F(B.x[j], B.v[j], t + h, C, A, B, j);
// Body C
m4[2][j] = h*(C.v[j]);
k4[2][j] = h*F(C.x[j], C.v[j], t + h, A, B, C, j);
}
for(int j=0; j<DIM-1; j++){
xA[j] += (m1[0][j] + 2*m2[0][j] + 2*m3[0][j] + m4[0][j])/6;
vA[j] += (k1[0][j] + 2*k2[0][j] + 2*k3[0][j] + k4[0][j])/6;
xB[j] += (m1[1][j] + 2*m2[1][j] + 2*m3[1][j] + m4[1][j])/6;
vB[j] += (k1[1][j] + 2*k2[1][j] + 2*k3[1][j] + k4[1][j])/6;
xC[j] += (m1[2][j] + 2*m2[2][j] + 2*m3[2][j] + m4[2][j])/6;
vC[j] += (k1[2][j] + 2*k2[2][j] + 2*k3[2][j] + k4[2][j])/6;
}
for(int j=0; j<DIM-1;j++){
A.v[j] = vA[j];
B.v[j] = vB[j];
C.v[j] = vC[j];
A.x[j] = xA[j];
B.x[j] = xB[j];
C.x[j] = xC[j];
}
#include <iostream>
#include <cmath>
#include <array>
#include <fstream>
#include "integrators.h"
#include <string>
#include <map>
// Value-Defintions of the different String values
enum StringValue { evNotDefined,
evStringValue1,
evStringValue2
};
// Map to associate the strings with the enum values
static std::map<std::string, StringValue> s_mapStringValues;
void Initialize(){
s_mapStringValues["euler"] = evStringValue1;
s_mapStringValues["rk4"] = evStringValue2;
}
static constexpr int DIM = 4;
static constexpr double G = 10;
static constexpr int N_BODIES = 3;
static constexpr int N_STEPS = 300000;
double distance(std::array<double, 3> r1, std::array<double, 3> r2){
return sqrt(pow(r1[0]-r2[0],2)+pow(r1[1]-r2[1],2)+pow(r1[2]-r2[2],2));
}
class Planet{
public:
double m;
std::array <double, 3> x;
std::array <double, 3> v;
std::array <double, 3> a;
double energy;
Planet (double mass, double x_position, double y_position, double z_position, double x_velocity, double y_velocity, double z_velocity) {
m = mass;
x[0] = x_position;
x[1] = y_position;
x[2] = z_position;
v[0] = x_velocity;
v[1] = y_velocity;
v[2] = z_velocity;
};
void computeKineticEnergy(){
energy = 0.5 * m * (pow(v[0],2)+pow(v[1],2)+pow(v[2],2));
}
};
double computePotentialEnergy(Planet planet1, Planet planet2, Planet planet3){
return -1 * G * ( planet3.m * planet1.m / distance(planet3.x, planet1.x) + planet3.m * planet2.m / distance(planet3.x, planet2.x) + planet1.m * planet2.m / distance(planet1.x, planet2.x));
}
double acceleration(Planet A, Planet B, Planet C, int axe){
// Compute the acceleration of the body C, specifying the axis.
double mass_A = A.m;
double mass_B = B.m;
double posx_A = A.x[0];
double posx_B = B.x[0];
double posx_C = C.x[0];
double posy_A = A.x[1];
double posy_B = B.x[1];
double posy_C = C.x[1];
double posz_A = A.x[2];
double posz_B = B.x[2];
double posz_C = C.x[2];
if (axe == 0){
return (-1 * G * (mass_A * (posx_C-posx_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posx_C-posx_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (axe == 1){
return (-1 * G * (mass_A * (posy_C-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posy_C-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (axe == 2){
return (-1 * G * (mass_A * (posz_C-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posz_C-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}
}
double F(double x, double v, double t, Planet A, Planet B, Planet C, int j ){
// Function to integrate via Runge-Kutta.
double mass_A = A.m;
double mass_B = B.m;
double posx_A = A.x[0];
double posx_B = B.x[0];
double posx_C = C.x[0];
double posy_A = A.x[1];
double posy_B = B.x[1];
double posy_C = C.x[1];
double posz_A = A.x[2];
double posz_B = B.x[2];
double posz_C = C.x[2];
if (j == 0){
return (-1 * G * (mass_A * (x-posx_A) / pow(sqrt(pow(x-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posx_B) / pow(sqrt(pow(x-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 1) {
return (-1 * G * (mass_A * (x-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(x-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(x-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 2) {
return (-1 * G * (mass_A * (x-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(x-posz_A,2)), 3) + mass_B * (x-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(x-posz_B,2)), 3)));
}
}
int main(int argc, char** argv){
double h = 0.0005;
Planet A(100, -10, 10, -11, -3, 0, 0);
Planet B(100, 0, 0, 0, 3, 0, 0);
Planet C(0, 10, 14, 12, 3, 0, 0);
double x_A[DIM][3];
double x_B[DIM][3];
double x_C[DIM][3];
double vx_A;
double vy_A;
double vz_A;
double vx_B;
double vy_B;
double vz_B;
double vx_C;
double vy_C;
double vz_C;
double mass_A = A.m;
double mass_B = B.m;
double mass_C = C.m;
x_A[0][0] = A.x[0];
x_B[0][0] = B.x[0];
x_C[0][0] = C.x[0];
vx_A = A.v[0];
vx_B = B.v[0];
vx_C = C.v[0];
x_A[1][0] = A.x[1];
x_B[1][0] = B.x[1];
x_C[1][0] = C.x[1];
vy_A = A.v[1];
vy_B = B.v[1];
vy_C = C.v[1];
x_A[2][0] = A.x[2];
x_B[2][0] = B.x[2];
x_C[2][0] = C.x[2];
vz_A = A.v[2];
vz_B = B.v[2];
vz_C = C.v[2];
Initialize();
std::ofstream file_energy("Total_energy_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_A("positions_A_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_B("positions_B_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_C("positions_C_" + std::string(argv[1]) + ".csv");
output_file_A<<"x;y;z"<<std::endl;
output_file_B<<"x;y;z"<<std::endl;
output_file_C<<"x;y;z"<<std::endl;
file_energy<<"k;p"<<std::endl;
double m1;
double k1;
double m2;
double k2;
double m3;
double k3;
double m4;
double k4;
std::array<double,3> vA = A.v; //condizione iniziale velocita
std::array<double,3> xA = A.x;
std::array<double,3> vB = B.v; //condizione iniziale velocita
std::array<double,3> xB = B.x;
std::array<double,3> vC = C.v; //condizione iniziale velocita
std::array<double,3> xC = C.x;
double t;
std::array<double,3> cm = computeCM(A, B, C);
if (argc>=2){
switch (s_mapStringValues[argv[1]]){
case evStringValue1:
// ==========================================================
// EULER
// ==========================================================
for (int i=0; i<N_STEPS-1; i++){
output_file_A << A.x[0] << ";" << A.x[1] << ";" << A.x[2]<< std::endl;
output_file_B << B.x[0] << ";" << B.x[1] << ";" << B.x[2]<< std::endl;
output_file_C << C.x[0] << ";" << C.x[1] << ";" << C.x[2]<< std::endl;
for(int j=0; j<DIM-1; j++){
A.a[j] = acceleration(B, C, A, j);
B.a[j] = acceleration(A, C, B, j);
C.a[j] = acceleration(B, A, C, j);
A.v[j] += A.a[j] * h;
B.v[j] += B.a[j] * h;
C.v[j] += C.a[j] * h;
x_A[j][0] = x_A[j][0] + A.v[j] * h;
x_B[j][0] = x_B[j][0] + B.v[j] * h;
x_C[j][0] = x_C[j][0] + C.v[j] * h;
}
for (int j=0;j<DIM-1;j++){
A.x[j] = x_A[j][0];
B.x[j] = x_B[j][0];
C.x[j] = x_C[j][0];
}
A.computeKineticEnergy();
B.computeKineticEnergy();
C.computeKineticEnergy();
cm = computeCM(A,B,C);
file_energy<<A.energy + B.energy + C.energy<<";"<< computePotentialEnergy(A, B, C)<<std::endl;
}
break;
case evStringValue2:{
// ==========================================================
// RUNGE KUTTA 4
// ==========================================================
for(int i=0; i<N_STEPS-1; i++){
for(int j=0; j<DIM-1; j++){
// body A
m1 = h*vA[j];
k1 = h*F(xA[j], vA[j], t, C, B, A, j);
m2 = h*(vA[j] + 0.5*k1);
k2 = h*F(xA[j]+0.5*m1, vA[j]+0.5*k1, t+0.5*h, C, B, A, j);
m3 = h*(vA[j] + 0.5*k2);
k3 = h*F(xA[j]+0.5*m2, vA[j]+0.5*k2, t+0.5*h, C, B, A, j);
m4 = h*(vA[j] + k3);
k4 = h*F(xA[j]+m3, vA[j]+k3, t+h, C, B, A, j);
xA[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vA[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1; j++){
// body B
m1 = h*vB[j];
k1 = h*F(xB[j], vB[j], t, A, C, B, j); //(x, v, t)
m2 = h*(vB[j] + 0.5*k1);
k2 = h*F(xB[j]+0.5*m1, vB[j]+0.5*k1, t+0.5*h, A, C, B, j);
m3 = h*(vB[j] + 0.5*k2);
k3 = h*F(xB[j]+0.5*m2, vB[j]+0.5*k2, t+0.5*h, A, C, B, j);
m4 = h*(vB[j] + k3);
k4 = h*F(xB[j]+m3, vB[j]+k3, t+h, A, C, B, j);
xB[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vB[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1; j++){
// body C
m1 = h*vC[j];
k1 = h*F(xC[j], vC[j], t, A, B, C, j); //(x, v, t)
m2 = h*(vC[j] + 0.5*k1);
k2 = h*F(xC[j]+0.5*m1, vC[j]+0.5*k1, t+0.5*h, A, B, C, j);
m3 = h*(vC[j] + 0.5*k2);
k3 = h*F(xC[j]+0.5*m2, vC[j]+0.5*k2, t+0.5*h, A, B, C, j);
m4 = h*(vC[j] + k3);
k4 = h*F(xC[j]+m3, vC[j]+k3, t+h, A, B, C, j);
xC[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vC[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1;j++){
A.v[j] = vA[j];
B.v[j] = vB[j];
C.v[j] = vC[j];
A.x[j] = xA[j];
B.x[j] = xB[j];
C.x[j] = xC[j];
}
output_file_A << A.x[0] << ";" << A.x[1] << ";" << A.x[2]<< std::endl;
output_file_B << B.x[0] << ";" << B.x[1] << ";" << B.x[2]<< std::endl;
output_file_C << C.x[0] << ";" << C.x[1] << ";" << C.x[2]<< std::endl;
A.computeKineticEnergy();
B.computeKineticEnergy();
C.computeKineticEnergy();
file_energy<<A.energy + B.energy + C.energy <<";"<< computePotentialEnergy(A, B, C)<<std::endl;
}
}
break;
default:
std::cout<<"Insert an argument between: euler, or rk4"<<std::endl;
return 0;
}
}else{
std::cout<<"Insert an argument between: euler, or rk4"<<std::endl;
return 0;
}
//----------------------------------------------------------------
output_file_A.close();
output_file_B.close();
output_file_C.close();
file_energy.close();
return 0;
}
#include <iostream>
#include <cmath>
#include <array>
#include <fstream>
#include "integrators.h"
#include <string>
#include <map>
// Value-Defintions of the different String values
enum StringValue { evNotDefined,
evStringValue1,
evStringValue2
};
// Map to associate the strings with the enum values
static std::map<std::string, StringValue> s_mapStringValues;
void Initialize(){
s_mapStringValues["euler"] = evStringValue1;
s_mapStringValues["rk4"] = evStringValue2;
}
static constexpr int DIM = 4;
static constexpr double G = 10;
static constexpr int N_BODIES = 3;
static constexpr int N_STEPS = 300000;
double distance(std::array<double, 3> r1, std::array<double, 3> r2){
return sqrt(pow(r1[0]-r2[0],2)+pow(r1[1]-r2[1],2)+pow(r1[2]-r2[2],2));
}
class Planet{
public:
double m;
std::array <double, 3> x;
std::array <double, 3> v;
std::array <double, 3> a;
double energy;
Planet (double mass, double x_position, double y_position, double z_position, double x_velocity, double y_velocity, double z_velocity) {
m = mass;
x[0] = x_position;
x[1] = y_position;
x[2] = z_position;
v[0] = x_velocity;
v[1] = y_velocity;
v[2] = z_velocity;
};
void computeKineticEnergy(){
energy = 0.5 * m * (pow(v[0],2)+pow(v[1],2)+pow(v[2],2));
}
};
double computePotentialEnergy(Planet planet1, Planet planet2, Planet planet3){
return -1 * G * ( planet3.m * planet1.m / distance(planet3.x, planet1.x) + planet3.m * planet2.m / distance(planet3.x, planet2.x) + planet1.m * planet2.m / distance(planet1.x, planet2.x));
}
double acceleration(Planet A, Planet B, Planet C, int axe){
// Compute the acceleration of the body C, specifying the axis.
double mass_A = A.m;
double mass_B = B.m;
double posx_A = A.x[0];
double posx_B = B.x[0];
double posx_C = C.x[0];
double posy_A = A.x[1];
double posy_B = B.x[1];
double posy_C = C.x[1];
double posz_A = A.x[2];
double posz_B = B.x[2];
double posz_C = C.x[2];
if (axe == 0){
return (-1 * G * (mass_A * (posx_C-posx_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posx_C-posx_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (axe == 1){
return (-1 * G * (mass_A * (posy_C-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posy_C-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (axe == 2){
return (-1 * G * (mass_A * (posz_C-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (posz_C-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}
}
double F(double x, double v, double t, Planet A, Planet B, Planet C, int j ){
// Function to integrate via Runge-Kutta.
double mass_A = A.m;
double mass_B = B.m;
double posx_A = A.x[0];
double posx_B = B.x[0];
double posx_C = C.x[0];
double posy_A = A.x[1];
double posy_B = B.x[1];
double posy_C = C.x[1];
double posz_A = A.x[2];
double posz_B = B.x[2];
double posz_C = C.x[2];
if (j == 0){
return (-1 * G * (mass_A * (x-posx_A) / pow(sqrt(pow(x-posx_A,2)+pow(posy_C-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posx_B) / pow(sqrt(pow(x-posx_B,2)+pow(posy_C-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 1) {
return (-1 * G * (mass_A * (x-posy_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(x-posy_A,2)+pow(posz_C-posz_A,2)), 3) + mass_B * (x-posy_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(x-posy_B,2)+pow(posz_C-posz_B,2)), 3)));
}else if (j == 2) {
return (-1 * G * (mass_A * (x-posz_A) / pow(sqrt(pow(posx_C-posx_A,2)+pow(posy_C-posy_A,2)+pow(x-posz_A,2)), 3) + mass_B * (x-posz_B) / pow(sqrt(pow(posx_C-posx_B,2)+pow(posy_C-posy_B,2)+pow(x-posz_B,2)), 3)));
}
}
int main(int argc, char** argv){
double h = 0.0005;
Planet A(100, -10, 10, -11, -3, 0, 0);
Planet B(100, 0, 0, 0, 3, 0, 0);
Planet C(0, 10, 14, 12, 3, 0, 0);
double x_A[DIM][3];
double x_B[DIM][3];
double x_C[DIM][3];
double vx_A;
double vy_A;
double vz_A;
double vx_B;
double vy_B;
double vz_B;
double vx_C;
double vy_C;
double vz_C;
double mass_A = A.m;
double mass_B = B.m;
double mass_C = C.m;
x_A[0][0] = A.x[0];
x_B[0][0] = B.x[0];
x_C[0][0] = C.x[0];
vx_A = A.v[0];
vx_B = B.v[0];
vx_C = C.v[0];
x_A[1][0] = A.x[1];
x_B[1][0] = B.x[1];
x_C[1][0] = C.x[1];
vy_A = A.v[1];
vy_B = B.v[1];
vy_C = C.v[1];
x_A[2][0] = A.x[2];
x_B[2][0] = B.x[2];
x_C[2][0] = C.x[2];
vz_A = A.v[2];
vz_B = B.v[2];
vz_C = C.v[2];
Initialize();
std::ofstream file_energy("Total_energy_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_A("positions_A_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_B("positions_B_" + std::string(argv[1]) + ".csv");
std::ofstream output_file_C("positions_C_" + std::string(argv[1]) + ".csv");
output_file_A<<"x;y;z"<<std::endl;
output_file_B<<"x;y;z"<<std::endl;
output_file_C<<"x;y;z"<<std::endl;
file_energy<<"k;p"<<std::endl;
double m1;
double k1;
double m2;
double k2;
double m3;
double k3;
double m4;
double k4;
std::array<double,3> vA = A.v; //condizione iniziale velocita
std::array<double,3> xA = A.x;
std::array<double,3> vB = B.v; //condizione iniziale velocita
std::array<double,3> xB = B.x;
std::array<double,3> vC = C.v; //condizione iniziale velocita
std::array<double,3> xC = C.x;
double t;
std::array<double,3> cm = computeCM(A, B, C);
if (argc>=2){
switch (s_mapStringValues[argv[1]]){
case evStringValue1:
// ==========================================================
// EULER
// ==========================================================
for (int i=0; i<N_STEPS-1; i++){
output_file_A << A.x[0] << ";" << A.x[1] << ";" << A.x[2]<< std::endl;
output_file_B << B.x[0] << ";" << B.x[1] << ";" << B.x[2]<< std::endl;
output_file_C << C.x[0] << ";" << C.x[1] << ";" << C.x[2]<< std::endl;
for(int j=0; j<DIM-1; j++){
A.a[j] = acceleration(B, C, A, j);
B.a[j] = acceleration(A, C, B, j);
C.a[j] = acceleration(B, A, C, j);
A.v[j] += A.a[j] * h;
B.v[j] += B.a[j] * h;
C.v[j] += C.a[j] * h;
x_A[j][0] = x_A[j][0] + A.v[j] * h;
x_B[j][0] = x_B[j][0] + B.v[j] * h;
x_C[j][0] = x_C[j][0] + C.v[j] * h;
}
for (int j=0;j<DIM-1;j++){
A.x[j] = x_A[j][0];
B.x[j] = x_B[j][0];
C.x[j] = x_C[j][0];
}
A.computeKineticEnergy();
B.computeKineticEnergy();
C.computeKineticEnergy();
cm = computeCM(A,B,C);
file_energy<<A.energy + B.energy + C.energy<<";"<< computePotentialEnergy(A, B, C)<<std::endl;
}
break;
case evStringValue2:{
// ==========================================================
// RUNGE KUTTA 4
// ==========================================================
for(int i=0; i<N_STEPS-1; i++){
for(int j=0; j<DIM-1; j++){
// body A
m1 = h*vA[j];
k1 = h*F(xA[j], vA[j], t, C, B, A, j);
m2 = h*(vA[j] + 0.5*k1);
k2 = h*F(xA[j]+0.5*m1, vA[j]+0.5*k1, t+0.5*h, C, B, A, j);
m3 = h*(vA[j] + 0.5*k2);
k3 = h*F(xA[j]+0.5*m2, vA[j]+0.5*k2, t+0.5*h, C, B, A, j);
m4 = h*(vA[j] + k3);
k4 = h*F(xA[j]+m3, vA[j]+k3, t+h, C, B, A, j);
xA[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vA[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1; j++){
// body B
m1 = h*vB[j];
k1 = h*F(xB[j], vB[j], t, A, C, B, j); //(x, v, t)
m2 = h*(vB[j] + 0.5*k1);
k2 = h*F(xB[j]+0.5*m1, vB[j]+0.5*k1, t+0.5*h, A, C, B, j);
m3 = h*(vB[j] + 0.5*k2);
k3 = h*F(xB[j]+0.5*m2, vB[j]+0.5*k2, t+0.5*h, A, C, B, j);
m4 = h*(vB[j] + k3);
k4 = h*F(xB[j]+m3, vB[j]+k3, t+h, A, C, B, j);
xB[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vB[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1; j++){
// body C
m1 = h*vC[j];
k1 = h*F(xC[j], vC[j], t, A, B, C, j); //(x, v, t)
m2 = h*(vC[j] + 0.5*k1);
k2 = h*F(xC[j]+0.5*m1, vC[j]+0.5*k1, t+0.5*h, A, B, C, j);
m3 = h*(vC[j] + 0.5*k2);
k3 = h*F(xC[j]+0.5*m2, vC[j]+0.5*k2, t+0.5*h, A, B, C, j);
m4 = h*(vC[j] + k3);
k4 = h*F(xC[j]+m3, vC[j]+k3, t+h, A, B, C, j);
xC[j] += (m1 + 2*m2 + 2*m3 + m4)/6;
vC[j] += (k1 + 2*k2 + 2*k3 + k4)/6;
}
for(int j=0; j<DIM-1;j++){
A.v[j] = vA[j];
B.v[j] = vB[j];
C.v[j] = vC[j];
A.x[j] = xA[j];
B.x[j] = xB[j];
C.x[j] = xC[j];
}
output_file_A << A.x[0] << ";" << A.x[1] << ";" << A.x[2]<< std::endl;
output_file_B << B.x[0] << ";" << B.x[1] << ";" << B.x[2]<< std::endl;
output_file_C << C.x[0] << ";" << C.x[1] << ";" << C.x[2]<< std::endl;
A.computeKineticEnergy();
B.computeKineticEnergy();
C.computeKineticEnergy();
file_energy<<A.energy + B.energy + C.energy <<";"<< computePotentialEnergy(A, B, C)<<std::endl;
}
}
break;
default:
std::cout<<"Insert an argument between: euler, or rk4"<<std::endl;
return 0;
}
}else{
std::cout<<"Insert an argument between: euler, or rk4"<<std::endl;
return 0;
}
//----------------------------------------------------------------
output_file_A.close();
output_file_B.close();
output_file_C.close();
file_energy.close();
return 0;
}
edit:
I think I have understood the problem. I have modified the code as follow, but now the energy decreases a lot and the trajectories become more and more little.
for(int i=0; i<N_STEPS-1; i++){
for(int j=0; j<DIM-1; j++){
// body A
m1[0][j] = h*vA[j];
k1[0][j] = h*F(xA[j], vA[j], t, C, B, A, j);
// body B
m1[1][j] = h*vA[j];
k1[1][j] = h*F(xA[j], vA[j], t, C, A, B, j);
// body C
m1[2][j] = h*vA[j];
k1[2][j] = h*F(xA[j], vA[j], t, A, B, C, j);
}
for(int j=0; j<DIM-1;j++){
A.v[j] += 0.5 * k1[0][j];
B.v[j] += 0.5 * k1[1][j];
C.v[j] += 0.5 * k1[2][j];
A.x[j] += 0.5 * m1[0][j];
B.x[j] += 0.5 * m1[1][j];
C.x[j] += 0.5 * m1[2][j];
}
for(int j=0; j<DIM-1; j++){
//Body A
m2[0][j] = h*(A.v[j]);
k2[0][j] = h*F(A.x[j], A.v[j], t+0.5*h, C, B, A, j);
//Body B
m2[1][j] = h*(B.v[j]);
k2[1][j] = h*F(B.x[j], B.v[j], t+0.5*h, C, A, B, j);
// Body C
m2[2][j] = h*(C.v[j]);
k2[2][j] = h*F(C.x[j], C.v[j], t+0.5*h, A, B, C, j);
}
for(int j=0; j<DIM-1;j++){
A.v[j] += 0.5 * k2[0][j];
B.v[j] += 0.5 * k2[1][j];
C.v[j] += 0.5 * k2[2][j];
A.x[j] += 0.5 * m2[0][j];
B.x[j] += 0.5 * m2[1][j];
C.x[j] += 0.5 * m2[2][j];
}
for(int j=0; j<DIM-1; j++){
//Body A
m3[0][j] = h*(A.v[j]);
k3[0][j] = h*F(A.x[j], A.v[j], t+0.5*h, C, B, A, j);
//Body B
m3[1][j] = h*(B.v[j]);
k3[1][j] = h*F(B.x[j], B.v[j], t+0.5*h, C, A, B, j);
// Body C
m3[2][j] = h*(C.v[j]);
k3[2][j] = h*F(C.x[j], C.v[j], t+0.5*h, A, B, C, j);
}
for(int j=0; j<DIM-1;j++){
A.v[j] += k3[0][j];
B.v[j] += k3[1][j];
C.v[j] += k3[2][j];
A.x[j] += m3[0][j];
B.x[j] += m3[1][j];
C.x[j] += m3[2][j];
}
for(int j=0; j<DIM-1; j++){
//Body A
m4[0][j] = h*(A.v[j]);
k4[0][j] = h*F(A.x[j], A.v[j], t + h, C, B, A, j);
//Body B
m4[1][j] = h*(B.v[j]);
k4[1][j] = h*F(B.x[j], B.v[j], t + h, C, A, B, j);
// Body C
m4[2][j] = h*(C.v[j]);
k4[2][j] = h*F(C.x[j], C.v[j], t + h, A, B, C, j);
}
for(int j=0; j<DIM-1; j++){
xA[j] += (m1[0][j] + 2*m2[0][j] + 2*m3[0][j] + m4[0][j])/6;
vA[j] += (k1[0][j] + 2*k2[0][j] + 2*k3[0][j] + k4[0][j])/6;
xB[j] += (m1[1][j] + 2*m2[1][j] + 2*m3[1][j] + m4[1][j])/6;
vB[j] += (k1[1][j] + 2*k2[1][j] + 2*k3[1][j] + k4[1][j])/6;
xC[j] += (m1[2][j] + 2*m2[2][j] + 2*m3[2][j] + m4[2][j])/6;
vC[j] += (k1[2][j] + 2*k2[2][j] + 2*k3[2][j] + k4[2][j])/6;
}
for(int j=0; j<DIM-1;j++){
A.v[j] = vA[j];
B.v[j] = vB[j];
C.v[j] = vC[j];
A.x[j] = xA[j];
B.x[j] = xB[j];
C.x[j] = xC[j];
}
Copy edited (e.g. ref. <https://en.wikipedia.org/wiki/Three-body_problem>, <https://en.wikipedia.org/wiki/Two-body_problem>, <https://en.wikipedia.org/wiki/C%2B%2B>, and <https://www.youtube.com/watch?v=1Dax90QyXgI&t=17m54s>).
Loading
Copy edited (e.g. ref. <https://en.wikipedia.org/wiki/Three-body_problem>, <https://en.wikipedia.org/wiki/Two-body_problem>, <https://en.wikipedia.org/wiki/C%2B%2B>, and <https://www.youtube.com/watch?v=1Dax90QyXgI&t=17m54s>).
Loading
Became Hot Network Question