C++三体问题模拟程序在macOS与Windows上表现差异求助

我和同事用C编写三体问题模拟程序，用于对比不同积分方案，当前采用欧拉法（Euler's method），通过终端调用g编译。但同一程序在macOS和Windows 10上运行结果存在差异：Windows平台中某一时刻天体的加速度会变为无穷大，而macOS平台无此问题，所有数据存储于三个.csv文件中，代码如下：

#include <iostream>
#include <cmath>
#include <array>
#include <fstream>

#define DIM 4
#define G  10//6.67408e-11
#define N_BODIES 3
#define N_STEPS 1000


float acceleration(float mass_1, float mass_2, float pos_1, float pos_2, float pos_3){ 
    //compute the acceleration along one axis of the body 3
    return -1 * G * (mass_1 * (pos_3-pos_1) / pow(abs(pos_3-pos_1), 3) + mass_2 * (pos_3-pos_2) / pow(abs(pos_3-pos_2), 3));
}

int main(){
    
    double mass_1 = 10, mass_2 = 20, mass_3 = 30;                       
    double x0_1[DIM-1] = {-10, 10, -11}, x0_2[DIM-1] = {0, 0, 0}, x0_3[DIM-1] = {10, 14, 12};                                  
    double v0_1[DIM-1] = {-3, 0, 0}, v0_2[DIM-1] = {0, 0, 0}, v0_3[DIM-1] = {3, 0, 0};                                 
    double h = 0.0005;

    double time[N_STEPS], x_1[DIM][N_STEPS], x_2[DIM][N_STEPS], x_3[DIM][N_STEPS], v_1[DIM][N_STEPS], v_2[DIM][N_STEPS], v_3[DIM][N_STEPS], a_1[DIM][N_STEPS], a_2[DIM][N_STEPS], a_3[DIM][N_STEPS];
          

    x_1[0][0] = x0_1[0];
    x_2[0][0] = x0_2[0];
    x_3[0][0] = x0_3[0];

    v_1[0][0] = v0_1[0];
    v_2[0][0] = v0_2[0];
    v_3[0][0] = v0_3[0];

    x_1[1][0] = x0_1[1];
    x_2[1][0] = x0_2[1];
    x_3[1][0] = x0_3[1];

    v_1[1][0] = v0_1[1];
    v_2[1][0] = v0_2[1];
    v_3[1][0] = v0_3[1];

    x_1[2][0] = x0_1[2];
    x_2[2][0] = x0_2[2];
    x_3[2][0] = x0_3[2];

    v_1[2][0] = v0_1[2];
    v_2[2][0] = v0_2[2];
    v_3[2][0] = v0_3[2];


//Function for the Euler method

for (int i=0; i<N_STEPS-1; i++){
    for(int j=0; j<DIM-1; j++){
    a_1[j][i] = acceleration(mass_2, mass_3, x_2[j][i], x_3[j][i], x_1[j][i]);
    a_2[j][i] = acceleration(mass_1, mass_3, x_1[j][i], x_3[j][i], x_2[j][i]);
    a_3[j][i] = acceleration(mass_1, mass_2, x_1[j][i], x_2[j][i], x_3[j][i]);
    
    v_1[j][i + 1] = v_1[j][i] + a_1[j][i] * h;
    v_2[j][i + 1] = v_2[j][i] + a_2[j][i] * h;
    v_3[j][i + 1] = v_3[j][i] + a_3[j][i] * h;
    
    x_1[j][i + 1] = x_1[j][i] + v_1[j][i] * h;
    x_2[j][i + 1] = x_2[j][i] + v_2[j][i] * h;
    x_3[j][i + 1] = x_3[j][i] + v_3[j][i] * h;
    
    }
    
}

    
    std::ofstream output_file_A("positions_A.csv");
    std::ofstream output_file_B("positions_B.csv");
    std::ofstream output_file_C("positions_C.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;
    

    for(int i = 0; i<N_STEPS-1; i++){
        output_file_A << a_1[0][i] << ";" << a_1[1][i] << ";" << a_1[2][i]<< std::endl;
        output_file_B << a_2[0][i] << ";" << a_2[1][i] << ";" << a_2[2][i]<< std::endl;
        output_file_C << a_3[0][i] << ";" << a_3[1][i] << ";" << a_3[2][i]<< std::endl;
    }    
    output_file_A.close();
    output_file_B.close(); 
    output_file_C.close();

    return 0;
}

问题根源分析 #

• 浮点数类型不匹配：acceleration函数的参数和返回值使用float类型，但主函数中所有变量都是double类型。跨平台时，float与double的精度、底层计算逻辑差异会导致数值偏差，Windows平台可能更早触发天体位置接近的临界情况，使分母趋近于0。
• 位置差的计算隐患：当两个天体在某一轴上的位置差极小时，pow(abs(pos_3-pos_1), 3)会趋近于0，直接导致加速度溢出为无穷大。不同平台的浮点数运算精度差异，会让Windows更早进入这个临界状态。
• 欧拉法的局限性：欧拉法是一阶积分方法，误差会随步数累积，跨平台的浮点计算差异会被不断放大，最终导致结果出现明显分歧。

解决方案 #

• 统一浮点数类型：将acceleration函数的参数和返回值全部改为double，避免类型转换带来的精度损失：

  double acceleration(double mass_1, double mass_2, double pos_1, double pos_2, double pos_3){ 
      //compute the acceleration along one axis of the body 3
      return -1 * G * (mass_1 * (pos_3-pos_1) / pow(abs(pos_3-pos_1), 3) + mass_2 * (pos_3-pos_2) / pow(abs(pos_3-pos_2), 3));
  }
  

• 添加距离阈值保护：在计算加速度前，判断位置差的绝对值是否小于一个极小值（比如1e-6），避免分母趋近于0：

  double acceleration(double mass_1, double mass_2, double pos_1, double pos_2, double pos_3){ 
      double diff1 = pos_3 - pos_1;
      double dist1 = abs(diff1);
      if(dist1 < 1e-6){
          dist1 = 1e-6; // 避免分母趋近于0
      }
      double diff2 = pos_3 - pos_2;
      double dist2 = abs(diff2);
      if(dist2 < 1e-6){
          dist2 = 1e-6;
      }
      return -1 * G * (mass_1 * diff1 / pow(dist1, 3) + mass_2 * diff2 / pow(dist2, 3));
  }
  

• 优化积分方法：替换为更高阶的积分方法，比如龙格-库塔法（RK4），减少误差累积，提升跨平台结果的一致性。
• 统一编译选项：跨平台编译时使用相同的优化等级和浮点标准，比如添加-std=c++17 -O2，尽量让浮点计算行为保持一致。

内容的提问来源于stack exchange，提问作者jack23456