c - 三体p‌r‌o‌b‌l‌e‌m(Burrau's pythagorean p‌r‌o‌b‌l‌e‌m.)与Runge-Kutta 4积分法(c语言)

Question

我正在尝试解决称为 Burrau 毕达哥拉斯问题的三体问题(初始条件 m1=3、m2=4、m3=5、x1 = 1、y1 = 3、x2 = -2、y2 = - 1, x3 = 1, y3 = -1) 用 Runge-Kutta 4 方法求解微分方程，但是当我使用 GNUPlot 绘制结果时(水平轴上的 x 位置，垂直轴上的 y 位置)我得到的都是直线，而不是 3 条困惑的轨迹。我也尝试过使用 Euler-Cromer 方法，但我仍然得到直线，所以问题一定不是 Runge-Kutta 4 方法。可能跟加速有关？ (f 函数)。不管是什么问题，我都想不通。这是代码:

#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

您需要一次对所有 12 个变量求积分，因为它们的方程是耦合的。通过分离积分步骤，您可以改变物理场并减少方法的阶数。典型的后果是引入一些漂移，即动量不守恒。

事实上，该方法仍然是 1 阶一致的，也就是说，如果您足够小地减小步长，那么您很可能会接近一些真实的物理行为。但是，如果步长太小，则积分会淹没在累积的浮点噪声中。