Problème à trois corps (problème de pythagoricien de Burrau) avec méthode d'intégration de Runge-Kutta 4 (langage c)

Question

J'essaie de résoudre le problème à trois corps connu sous le nom de problème de pythagoricien de Burrau (conditions initiales m1 = 3, m2 = 4 , m3 = 5, x1 = 1, y1 = 3, x2 = -2, y2 = -1, x3 = 1, y3 = -1) avec la méthode de Runge-Kutta 4 pour résoudre les équations différentielles, mais quand je trace le les résultats utilisant GNUPlot (x positions sur l'axe horizontal, y positions sur l'axe vertical) tout ce que je reçois sont des lignes droites, au lieu de 3 trajectoires chaotiques comme il se doit. J'ai aussi essayé d'utiliser la méthode d'Euler-Cromer, et j'ai toujours des lignes droites, donc le problème ne doit pas être avec la méthode Runge-Kutta 4. Peut-être que c'est avec l'accélération? (f fonction). Quel que soit le problème, je ne peux pas le comprendre. Voici le 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; 
} 

Answer 1

Vous devez intégrer les 12 variables à la fois, car leurs équations sont couplées. En séparant les étapes d'intégration, vous changez la physique et vous réduisez l'ordre de la méthode. Une conséquence typique est l'introduction d'une certaine dérive, c'est-à-dire que l'élan n'est pas conservé. Comme c'est le cas, la méthode est toujours d'ordre 1, c'est-à-dire que si vous réduisez suffisamment la taille du pas, vous vous rapprocherez vraisemblablement d'un comportement physique réaliste. Cependant, si la taille du pas devient trop petite, l'intégration se noie dans le bruit de point flottant accumulé.