I'm trying to solve the Three-body problem known as Burrau's pythagorean problem (initial conditions m1=3, m2=4, m3=5, x1 = 1, y1 = 3, x2 = -2, y2 = -1, x3 = 1, y3 = -1) with the Runge-Kutta 4 method to solve the differential equations, but when I plot the results using GNUPlot (x positions on the horizontal axis, y positions on vertical axis) all I get are straight lines, instead of 3 chaotic trajectories as it should be. I've tried using the Euler-Cromer method too, and I still get straight lines, so the issue must not be with the Runge-Kutta 4 method. Maybe it's with the acceleration? (f function). Whatever the issue is, I can't figure it out. Here's the code:
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
/* N.B the units have been normalized so that the gravitational coefficient is 1 */
typedef struct {
long long unsigned int N;
double deltat;
double tmax;
double tn;
} variabile;
typedef struct {
double m;
double xn;
double yn;
double vxn;
double vyn;
} corpo;
double f(double, double, double, double, double, double, double, double);
void RungeKutta4(variabile*, corpo*, corpo*, corpo*, FILE*);
void Euler(variabile*, corpo*, corpo*, corpo*, FILE*);
int main(int argc, char ** argv) {
variabile variab;
variabile *var = &variab;
corpo body1, body2, body3;
corpo *c1 = &body1, *c2 = &body2, *c3 = &body3;
FILE *fp;
fp = fopen("tbp_RungeKutta4.dat", "w");
if (argc != 3) {
system("clear");
fprintf(stderr, "\nCorrect syntax: \n\n'tmax deltat'\n\n\n\n\n");
exit(EXIT_FAILURE);
}
variab.tmax = atof(argv[1]);
variab.deltat = atof(argv[2]);
variab.N = (variab.tmax)/(variab.deltat);
/* Initial conditions for body1: */
body1.m = 3.0;
body1.xn = 1.0;
body1.yn = 3.0;
body1.vxn = 0.0;
body1.vyn = 0.0;
/* Initial conditions for body2: */
body2.m = 4.0;
body2.xn = -2.0;
body2.yn = -1.0;
body2.vxn = 0.0;
body2.vyn = 0.0;
/* Initial conditions for body3: */
body3.m = 5.0;
body3.xn = 1.0;
body3.yn = -1.0;
body3.vxn = 0.0;
body3.vyn = 0.0;
system("clear");
printf("Calculating...\n");
RungeKutta4(var, c1, c2, c3, fp);
printf("...Success!\n");
return EXIT_SUCCESS;
} /* END OF MAIN */
double f(double xi, double yi, double mj, double xj, double yj, double mk, double xk, double yk) {
double a;
a = (xj-xi)*(mj/pow((xi-xj)*(xi-xj) + (yi-yj)*(yi-yj), 1.5)) + (xk-xi)*(mk/pow((xi-xk)*(xi-xk) + (yi-yk)*(yi-yk), 1.5));
return a;
}
void RungeKutta4(variabile *var, corpo *c1, corpo *c2, corpo *c3, FILE *fp) {
int i;
long long unsigned int N = var->N;
double tn = var->tn, dt = var->deltat;
double m1, xn1, yn1, vxn1, vyn1;
double m2, xn2, yn2, vxn2, vyn2;
double m3, xn3, yn3, vxn3, vyn3;
double dp1, dp2, dp3, dp4, dv1, dv2, dv3, dv4;
double temp_dxn1, temp_dvxn1, temp_dyn1, temp_dvyn1;
double temp_dxn2, temp_dvxn2, temp_dyn2, temp_dvyn2;
double temp_dxn3, temp_dvxn3, temp_dyn3, temp_dvyn3;
/* Temporary variables for body1 */
m1 = c1->m;
xn1 = c1->xn;
yn1 = c1->yn;
vxn1 = c1->vxn;
vyn1 = c1->vyn;
/* Temporary variables for body2 */
m2 = c2->m;
xn2 = c2->xn;
yn2 = c2->yn;
vxn2 = c2->vxn;
vyn2 = c2->vyn;
/* Temporary variables for body3 */
m3 = c3->m;
xn3 = c3->xn;
yn3 = c3->yn;
vxn3 = c3->vxn;
vyn3 = c3->vyn;
/* Writing the initial values (time=0) on file */
fprintf(fp, "#x1\t\ty1\t\tx2\t\ty2\t\tx3\t\ty3\t\ttempo\n");
fprintf(fp, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", xn1, yn1, xn2, yn2, xn3, yn3, tn);
for (i=1; i<=N; i++) {
tn += dt;
/* BODY c1: */
dp1 = vxn1*dt;
dv1 = f(xn1, yn1, m2, xn2, yn2, m3, xn3, yn3)*dt;
dp2 = (vxn1 + 0.5*dv1)*dt;
dv2 = f(xn1 + 0.5*dp1, yn1, m2, xn2, yn2, m3, xn3, yn3)*dt;
dp3 = (vxn1 + 0.5*dv2)*dt;
dv3 = f(xn1 + 0.5*dp2, yn1, m2, xn2, yn2, m3, xn3, yn3)*dt;
dp4 = (vxn1 + dv3)*dt;
dv4 = f(xn1 + dp3, yn1, m2, xn2, yn2, m3, xn3, yn3)*dt;
temp_dxn1 = (1./6.)*(dp1 + 2*dp2 + 2*dp3 + dp4);
temp_dvxn1 = (1./6.)*(dv1 + 2*dv2 + 2*dv3 + dv4);
dp1 = vyn1*dt;
dv1 = f(yn1, xn1, m2, yn2, xn2, m3, yn3, xn3)*dt;
dp2 = (vyn1 + 0.5*dv1)*dt;
dv2 = f(yn1 + 0.5*dp1, xn1, m2, yn2, xn2, m3, yn3, xn3)*dt;
dp3 = (vyn1 + 0.5*dv2)*dt;
dv3 = f(yn1 + 0.5*dp2, xn1, m2, yn2, xn2, m3, yn3, xn3)*dt;
dp4 = (vyn1 + dv3)*dt;
dv4 = f(yn1 + dp3, xn1, m2, yn2, xn2, m3, yn3, xn3)*dt;
temp_dyn1 = (1./6.)*(dp1 + 2*dp2 + 2*dp3 + dp4);
temp_dvyn1 = (1./6.)*(dv1 + 2*dv2 + 2*dv3 + dv4);
/* BODY c2: */
dp1 = vxn2*dt;
dv1 = f(xn2, yn2, m1, xn1, yn1, m3, xn3, yn3)*dt;
dp2 = (vxn2 + 0.5*dv1)*dt;
dv2 = f(xn2 + 0.5*dp1, yn2, m1, xn1, yn1, m3, xn3, yn3)*dt;
dp3 = (vxn2 + 0.5*dv2)*dt;
dv3 = f(xn2 + 0.5*dp2, yn2, m1, xn1, yn1, m3, xn3, yn3)*dt;
dp4 = (vxn2 + dv3)*dt;
dv4 = f(xn2 + dp3, yn2, m1, xn1, yn1, m3, xn3, yn3)*dt;
temp_dxn2 = (1./6.)*(dp1 + 2*dp2 + 2*dp3 + dp4);
temp_dvxn2 = (1./6.)*(dv1 + 2*dv2 + 2*dv3 + dv4);
dp1 = vyn2*dt;
dv1 = f(yn2, xn2, m1, yn1, xn1, m3, yn3, xn3)*dt;
dp2 = (vyn2 + 0.5*dv1)*dt;
dv2 = f(yn2 + 0.5*dp1, xn2, m1, yn1, xn1, m3, yn3, xn3)*dt;
dp3 = (vyn2 + 0.5*dv2)*dt;
dv3 = f(yn2 + 0.5*dp2, xn2, m1, yn1, xn1, m3, yn3, xn3)*dt;
dp4 = (vyn2 + dv3)*dt;
dv4 = f(yn2 + dp3, xn2, m1, yn1, xn1, m3, yn3, xn3)*dt;
temp_dyn2 = (1./6.)*(dp1 + 2*dp2 + 2*dp3 + dp4);
temp_dvyn2 = (1./6.)*(dv1 + 2*dv2 + 2*dv3 + dv4);
/* BODY c3: */
dp1 = vxn3*dt;
dv1 = f(xn3, yn3, m1, xn1, yn1, m2, xn2, yn2)*dt;
dp2 = (vxn3 + 0.5*dv1)*dt;
dv2 = f(xn3 + 0.5*dp1, yn3, m1, xn1, yn1, m2, xn2, yn2)*dt;
dp3 = (vxn3 + 0.5*dv2)*dt;
dv3 = f(xn3 + 0.5*dp2, yn3, m1, xn1, yn1, m2, xn2, yn2)*dt;
dp4 = (vxn3 + dv3)*dt;
dv4 = f(xn3 + dp3, yn3, m1, xn1, yn1, m2, xn2, yn2)*dt;
temp_dxn3 = (1./6.)*(dp1 + 2*dp2 + 2*dp3 + dp4);
temp_dvxn3 = (1./6.)*(dv1 + 2*dv2 + 2*dv3 + dv4);
dp1 = vyn3*dt;
dv1 = f(yn3, xn3, m1, yn1, xn1, m2, yn2, xn2)*dt;
dp2 = (vyn3 + 0.5*dv1)*dt;
dv2 = f(yn3 + 0.5*dp1, xn3, m1, yn1, xn1, m2, yn2, xn2)*dt;
dp3 = (vyn3 + 0.5*dv2)*dt;
dv3 = f(yn3 + 0.5*dp2, xn3, m1, yn1, xn1, m2, yn2, xn2)*dt;
dp4 = (vyn3 + dv3)*dt;
dv4 = f(yn3 + dp3, xn3, m1, yn1, xn1, m2, yn2, xn2)*dt;
temp_dyn3 = (1./6.)*(dp1 + 2*dp2 + 2*dp3 + dp4);
temp_dvyn3 = (1./6.)*(dv1 + 2*dv2 + 2*dv3 + dv4);
/************** END OF RK4 **************/
/* Increasing variables */
xn1 += temp_dxn1;
vxn1 += temp_dvxn1;
yn1 += temp_dyn1;
vyn1 += temp_dvyn1;
xn2 += temp_dxn2;
vxn2 += temp_dvxn2;
yn2 += temp_dyn2;
vyn2 += temp_dvyn2;
xn3 += temp_dxn3;
vxn3 += temp_dvxn3;
yn3 += temp_dyn3;
vyn3 += temp_dvyn3;
/* Updating variables on the structs from the temporary variables */
var->tn = tn;
/* Body 1 */
c1->xn = xn1;
c1->yn = yn1;
c1->vxn = vxn1;
c1->vyn = vyn1;
/* Body 2 */
c2->xn = xn2;
c2->yn = yn2;
c2->vxn = vxn2;
c2->vyn = vyn2;
/* Body 3 */
c3->xn = xn3;
c3->yn = yn3;
c3->vxn = vxn3;
c3->vyn = vyn3;
/* Writing data on file */
fprintf(fp, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", xn1, yn1, xn2, yn2, xn3, yn3, tn);
}
return;
}
You need to integrate all 12 variables at once, since their equations are coupled. By separating the integration steps, you change the physics and you reduce the order of the method. A typical consequence is the introduction of some drift, i.e., the momentum is not conserved.
As it is, the method is still order 1 consistent, that is, if you reduce the step size enough, then in all probability you will come close to some realistic physical behaviour. However, if the step size gets too small, then the integration drowns in the accumulating floating point noise.