1

I'm trying to write a code to solve the n-body problem using Runge Kutta 4 integration algorithm. I'm testing the code using two bodies with masses uniformly distributed between 0 and 1, position distributed following a density law proportional to 1/r^2 and velocity distributed as Maxwell-Boltzmann distribution. I've tried to integrate the system for different tmax, but I get the orbit in the plots and I can't figure out the problem. Any help would be appreciated.

Here's the code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct {
    double x, y, z;
}vector;

double *masse;
vector *pos, *vel, *forza;



int main(int argc, char** argv){
    double epsilon, dt, tmax,t;
    double dx, dy, dz, dist, invdist, invdist3;
    int N, i, l, m, n, s, j;
    double a, b;
    vector *k1, *k2, *k3, *k4, *w1, *w2, *w3, *w4, *pos1, *vel1;
    if(argc!=5) {
        fprintf(stdout,"Il programma prende in input il softening, il passo d'integrazione, il tempo massimo d'integrazione e il numero di corpi del sistema\n", argv[0]);
        exit(1);
    }

    epsilon=strtod(argv[1],NULL);
    dt=strtod(argv[2],NULL);
    tmax=strtod(argv[3],NULL);
    N=strtod(argv[4],NULL);

    FILE* fp=fopen("Cond_ini.out", "r");
    if(fp==NULL){
        perror("Errore: file non trovato\n");
        exit(1);
    }

    masse=(double*)malloc(N*sizeof(double));
    pos=(vector*)malloc(N*sizeof(vector));
    vel=(vector*)malloc(N*sizeof(vector));
    forza=(vector*)malloc(N*sizeof(vector));
    k1=(vector*)malloc(N*sizeof(vector));
    k2=(vector*)malloc(N*sizeof(vector));
    k3=(vector*)malloc(N*sizeof(vector));
    k4=(vector*)malloc(N*sizeof(vector));
    w1=(vector*)malloc(N*sizeof(vector));
    w2=(vector*)malloc(N*sizeof(vector));
    w3=(vector*)malloc(N*sizeof(vector));
    w4=(vector*)malloc(N*sizeof(vector));
    pos1=(vector*)malloc(N*sizeof(vector));
    vel1=(vector*)malloc(N*sizeof(vector));





    for(i=0;i<N;i++){
        fscanf(fp,"%lf %lf %lf %lf %lf %lf %lf", &masse[i], &pos[i].x, &pos[i].y, &pos[i].z, &vel[i].x, &vel[i].y, &vel[i].z);
    }
    fclose(fp);

    printf("Condizioni iniziali:\n");
    for(l=0;l<N;l++){
        printf("%lf %lf %lf %lf %lf %lf %lf\n", masse[l], pos[l].x, pos[l].y, pos[l].z, vel[l].x, vel[l].y, vel[l].z);
    }

    for(t=0;t<tmax;t+=dt){
        for(m=0;m<N;m++){
            for(n=0;n<N;n++){
                if(m!=n){

                    k1[n].x=dt*vel[n].x;
                    k1[n].y=dt*vel[n].y;
                    k1[n].z=dt*vel[n].z;

                    dx=pos[n].x-pos[m].x;
                    dy=pos[n].y-pos[m].y;
                    dz=pos[n].z-pos[m].z;
                    dist=(dx*dx)+(dy*dy)+(dz*dz)+(epsilon*epsilon);
                    invdist=(1/sqrt(dist));
                    invdist3=(invdist*invdist*invdist);
                    forza[n].x=dx*invdist3*masse[n];
                    forza[n].y=dy*invdist3*masse[n];
                    forza[n].z=dz*invdist3*masse[n];

                    w1[n].x=dt*forza[n].x;
                    w1[n].y=dt*forza[n].y;
                    w1[n].z=dt*forza[n].z;

                    pos1[n].x=pos[n].x+(0.5*k1[n].x);
                    pos1[n].y=pos[n].y+(0.5*k1[n].y);
                    pos1[n].z=pos[n].z+(0.5*k1[n].z);
                    vel1[n].x=vel[n].x+(0.5*w1[n].x);
                    vel1[n].y=vel[n].y+(0.5*w1[n].y);
                    vel1[n].z=vel[n].z+(0.5*w1[n].z);

                    k2[n].x=dt*(vel[n].x+(0.5*w1[n].x));
                    k2[n].y=dt*(vel[n].y+(0.5*w1[n].y));
                    k2[n].z=dt*(vel[n].z+(0.5*w1[n].z));

                    dx=pos1[n].x-pos[m].x;
                    dy=pos1[n].y-pos[m].y;
                    dz=pos1[n].z-pos[m].z;
                    dist=(dx*dx)+(dy*dy)+(dz*dz)+(epsilon*epsilon);
                    invdist=(1/sqrt(dist));
                    invdist3=(invdist*invdist*invdist);
                    forza[n].x=dx*invdist3*masse[n];
                    forza[n].y=dy*invdist3*masse[n];
                    forza[n].z=dz*invdist3*masse[n];

                    w2[n].x=dt*forza[n].x;
                    w2[n].y=dt*forza[n].y;
                    w2[n].z=dt*forza[n].z;

                    pos1[n].x=pos[n].x+(0.5*k2[n].x);
                    pos1[n].y=pos[n].y+(0.5*k2[n].y);
                    pos1[n].z=pos[n].z+(0.5*k2[n].z);
                    vel1[n].x=vel[n].x+(0.5*w2[n].x);
                    vel1[n].y=vel[n].y+(0.5*w2[n].y);
                    vel1[n].z=vel[n].z+(0.5*w2[n].z);

                    k3[n].x=dt*(vel[n].x+(0.5*w2[n].x));
                    k3[n].y=dt*(vel[n].y+(0.5*w2[n].y));
                    k3[n].z=dt*(vel[n].z+(0.5*w2[n].z));

                    dx=pos1[n].x-pos[m].x;
                    dy=pos1[n].y-pos[m].y;
                    dz=pos1[n].z-pos[m].z;
                    dist=(dx*dx)+(dy*dy)+(dz*dz)+(epsilon*epsilon);
                    invdist=(1/sqrt(dist));
                    invdist3=(invdist*invdist*invdist);
                    forza[n].x=dx*invdist3*masse[n];
                    forza[n].y=dy*invdist3*masse[n];
                    forza[n].z=dy*invdist3*masse[n];

                    w3[n].x=dt*forza[n].x;
                    w3[n].y=dt*forza[n].y;
                    w3[n].z=dt*forza[n].z;

                    pos1[n].x=pos[n].x+(k3[n].x);
                    pos1[n].y=pos[n].y+(k3[n].y);
                    pos1[n].z=pos[n].z+(k3[n].z);
                    vel1[n].x=vel[n].x+(w3[n].x);
                    vel1[n].y=vel[n].y+(w3[n].y);
                    vel1[n].z=vel[n].z+(w3[n].z);


                    k4[n].x=dt*(vel[n].x+w3[n].x);
                    k4[n].y=dt*(vel[n].y+w3[n].y);
                    k4[n].z=dt*(vel[n].z+w3[n].z);

                    dx=pos1[n].x-pos[m].x;
                    dy=pos1[n].y-pos[m].y;
                    dz=pos1[n].z-pos[m].z;
                    dist=(dx*dx)+(dy*dy)+(dz*dz)+(epsilon*epsilon);
                    invdist=(1/sqrt(dist));
                    invdist3=(invdist*invdist*invdist);
                    forza[n].x=dx*invdist3*masse[n];
                    forza[n].y=dy*invdist3*masse[n];
                    forza[n].z=dy*invdist3*masse[n];

                    w4[n].x=dt*forza[n].x;
                    w4[n].y=dt*forza[n].y;
                    w4[n].z=dt*forza[n].z;

                    a=k1[n].x+(2*k2[n].x)+(2*k3[n].x)+k4[n].x;
                    a=a/6;
                    pos1[n].x=pos[n].x+a;
                    a=k1[n].y+(2*k2[n].y)+(2*k3[n].y)+k4[n].y;
                    a=a/6;
                    pos1[n].y=pos[n].y+a;
                    a=k1[n].z+(2*k2[n].z)+(2*k3[n].z)+k4[n].z;
                    a=a/6;
                    pos1[n].z=pos[n].z+a;
                    b=w1[n].x+(2*w2[n].x)+(2*w3[n].x)+w4[n].x;
                    b=b/6;
                    vel1[n].x=vel[n].x+a;
                    b=w1[n].y+(2*w2[n].y)+(2*w3[n].y)+w4[n].y;
                    b=b/6;
                    vel1[n].y=vel[n].y+a;
                    b=w1[n].z+(2*w2[n].z)+(2*w3[n].z)+w4[n].z;
                    b=b/6;
                    vel1[n].z=vel[n].z+a;
                }
            }
            for(j=0;j<N;j++) {
                forza[j].x=0;
                forza[j].y=0;
                forza[j].z=0;
            }
        }
        for(i=0;i<N;i++){
            pos[i].x=pos1[i].x;
            pos[i].y=pos1[i].y;
            pos[i].z=pos1[i].z;
            vel[i].x=vel1[i].x;
            vel[i].y=vel1[i].y;
            vel[i].z=vel1[i].z;
            printf("%lf %lf %lf %lf %lf %lf %lf\n", t, pos[i].x, pos[i].y, pos[i].z, vel[i].x, vel[i].y, vel[i].z);
            /*forza[i].x=0;
            forza[i].y=0;
            forza[i].z=0;*/
        }
    }
}

And here's a plot of orbits at different tmax:

tmax=5 tmax=2 tmax=3

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
Maere
  • 11
  • 1
  • I have tried to understand the problem you are facing, Is is a **C Programming** probelm or physics problem, what is **forza** – Kuldip Chaudhari Aug 23 '20 at 16:48
  • @KuldipChaudhari it's a physics problem related to the motions of objects that interact each other gravitationally – sebastian Aug 23 '20 at 17:01
  • @KuldipChaudhari "forza" is body's acceleration – Maere Aug 23 '20 at 17:13
  • Thanks, @Maere, I will try to solve this, I need time to understand this – Kuldip Chaudhari Aug 23 '20 at 17:17
  • regarding: `fprintf(stdout,"Il programma prende in input il softening, il passo d'integrazione, il tempo massimo d'integrazione e il numero di corpi del sistema\n", argv[0]);` this will fail to properly compile due to the parameter: `argv[0]` but no `%s` in the format string. suggest: `fprintf(stdout,"USAGE: %s Il programma prende in input il softening, il passo d'integrazione, il tempo massimo d'integrazione e il numero di corpi del sistema\n", argv[0]);` *untitled1.c:21:24: warning: too many arguments for format [-Wformat-extra-args]* – user3629249 Aug 24 '20 at 22:48
  • running the posted code through the compiler results in many warning messages `untitled1.c:41:25: warning: conversion to ‘long unsigned int’ from ‘int’ may change the sign of the result [-Wsign-conversion]` – user3629249 Aug 24 '20 at 22:51
  • OT: regarding statements like;: `masse=(double*)malloc(N*sizeof(double));` 1) In c, the returned type is `void*` which can be assigned to any pointer. Casting just clutters the code and is error prone. Suggest removing the cast. 2) always check (!=NULL) the returned value to assure the operation was successful. If not successful (==NULL) then call `perror( "malloc failed" );` to output both the error message and the text reason the system thinks the error occurred. – user3629249 Aug 24 '20 at 22:57
  • the posted code fails to pass all those allocated memory pointers to `free()` before exiting the program. The OS might cleanup for you, but this sloppiness will burn you when working with embedded systems. – user3629249 Aug 24 '20 at 23:01
  • OT: for ease of readability and understanding: 1) please follow the axiom: *only one statement per line and (at most) one variable declaration per statement.* 2) please use meaningful variable names. Names like: `vector *k1, *k2, *k3, *k4, *w1, *w2, *w3, *w4, *pos1, *vel1;` and `int N, i, l, m, n, s, j;`, etc are meaningless to us. – user3629249 Aug 24 '20 at 23:07
  • there is no/zero/nada reason to allocate all those dynamic memory allocations. Rather use the Variable Length Array feature of C. I.E. use something like: `double massie[ N * sizeof( vector ) ];` after having obtained `N` – user3629249 Aug 24 '20 at 23:09

1 Answers1

0

In looking through your code, your first error is an algorithmical one. The order of loops has to be time steps, then the RK4 stages, and inside that the summation of the interaction terms/forces of the stellar bodies. This means especially that you need to compute the interactions separately and anew for each stage of the Runge-Kutta method.

It could make sense to make a separate procedure that fills the force array from the position array to avoid copy-paste-edit errors.

For a well-documented test problem of this type see the Pleiades system of the IVP test suite https://archimede.dm.uniba.it/~testset/testsetivpsolvers/?page_id=26

Cf. also How can I update c++ class members with a function? where I posted code implementing these principles/corrections on this test problem for a simulation with the Verlet method (and in C++).

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51