如何解决三体问题模拟程序返回'inf'的问题？
2026-6-23
解决三体问题模拟中出现inf/nan的问题
我正在研究经典物理中的三体问题，编写了如下C++程序：
#include <math.h> #include <iostream> using namespace std; class System_variables { private: // 系统变量：位置、速度、加速度 long double body_1_pos[3]; long double body_2_pos[3]; long double body_3_pos[3]; long double body_1_v[3]; long double body_2_v[3]; long double body_3_v[3]; long double acc_body_1[3]; long double acc_body_2[3]; long double acc_body_3[3]; // 所需常量 const long double body_mass = 6 * pow(10, 24); double time_step; const long double G = 6.67 * pow(10, -11); public: System_variables(long double body_1_pos_i[3], long double body_2_pos_i[3], long double body_3_pos_i[3], long double body_1_v_i[3], long double body_2_v_i[3], long double body_3_v_i[3]) { // 初始化位置和速度 for (int m = 0; m<=2; m++) { body_1_pos[m] = body_1_pos_i[m]; body_2_pos[m] = body_2_pos_i[m]; body_3_pos[m] = body_3_pos_i[m]; body_1_v[m] = body_1_v_i[m]; body_2_v[m] = body_2_v_i[m]; body_3_v[m] = body_3_v_i[m]; } } void init_time_step(double step) { // 初始化时间步长 time_step = step; } void calculate_acc() { // 用牛顿方程计算天体加速度 for (int m=0; m<=2; m++) { acc_body_1[m] = G * body_mass/pow((body_2_pos[m] - body_1_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_1_pos[m]), 2); acc_body_2[m] = G * body_mass/pow((body_1_pos[m] - body_2_pos[m]), 2) + G * body_mass/pow((body_3_pos[m] - body_2_pos[m]), 2); acc_body_3[m] = G * body_mass/pow((body_1_pos[m] - body_3_pos[m]), 2) + G * body_mass/pow((body_2_pos[m] - body_3_pos[m]), 2); }; return; } void apply_acc_to_v() { // 将加速度应用到速度 for (int m=0; m<=2; m++) { body_1_v[m] += acc_body_1[m] * time_step; body_2_v[m] += acc_body_2[m] * time_step; body_3_v[m] += acc_body_3[m] * time_step; } return; } void apply_v_to_pos() { // 将速度应用到位置 for (int m=0; m<=2; m++) { body_1_pos[m] += body_1_v[m] * time_step; body_2_pos[m] += body_2_v[m] * time_step; body_3_pos[m] += body_3_v[m] * time_step; } return; } long double get_body_1_pos(int index) { // 返回天体1对应坐标轴的位置 return body_1_pos[index]; } long double get_body_2_pos(int index) { // 返回天体2对应坐标轴的位置 return body_2_pos[index]; } long double get_body_3_pos(int index) { // 返回天体3对应坐标轴的位置 return body_3_pos[index]; } long double get_body_1_v(int index) { // 返回天体1对应坐标轴的速度 return body_1_v[index]; } long double get_body_2_v(int index) { // 返回天体2对应坐标轴的速度 return body_2_v[index]; } long double get_body_3_v(int index) { // 返回天体3对应坐标轴的速度 return body_3_v[index]; } }; int main() { long double array_1[3] = {100,100,100}; long double array_2[3] ={10000000000, 10000000000, 10000000000}; long double array_3[3] = {10000000000, 10000000000, -10000000000}; long double array_0[3] = {0,0,0}; System_variables system(array_1, array_2, array_3, array_0, array_0, array_0); system.init_time_step(100); for (int time = 0; time <= 10; time++) { cout << "x position of body 1: " << system.get_body_1_pos(0) << endl << "y position of body 1: " << system.get_body_1_pos(1) << endl << "z position of body 1: " << system.get_body_1_pos(2) << endl << endl; cout << "x position of body 2: " << system.get_body_2_pos(0) << endl << "y position of body 2: " << system.get_body_2_pos(1) << endl << "z position of body 2: " << system.get_body_2_pos(2) << endl << endl; cout << "x position of body 3: " << system.get_body_3_pos(0) << endl << "y position of body 3: " << system.get_body_3_pos(1) << endl << "z position of body 3: " << system.get_body_3_pos(2) << endl << endl; system.calculate_acc(); system.apply_acc_to_v(); system.apply_v_to_pos(); } }
我尝试过用long double替代double，以为更高精度能解决运行后输出inf的问题，但没用；换成float也一样，输出结果如下：
x position of body 1: 100 y position of body 1: 100 z position of body 1: 100 x position of body 2: 1e+10 y position of body 2: 1e+10 z position of body 2: 1e+10 x position of body 3: 1e+10 y position of body 3: 1e+10 z position of body 3: -1e+10 x position of body 1: 100.08 y position of body 1: 100.08 z position of body 1: 100.08 x position of body 2: inf y position of body 2: inf z position of body 2: 1e+10 x position of body 3: inf y position of body 3: inf z position of body 3: -1e+10 x position of body 1: 100.16 y position of body 1: 100.16 z position of body 1: 100.24 x position of body 2: -nan(ind) y position of body 2: -nan(ind) z position of body 2: 1e+10 x position of body 3: -nan(ind) y position of body 3: -nan(ind) z position of body 3: -1e+10 x position of body 1: -nan(ind) y position of body 1: -nan(ind) z position of body 1: 100.48 x position of body 2: -nan(ind) y position of body 2: -nan(ind) z position of body 2: 1e+10 x position of body 3: -nan(ind) y position of body 3: -nan(ind) z position of body 3: -1e+10
我考虑过用指针，但不确定是不是正确方向，请问怎么解决inf的问题？
问题根源分析
- 引力加速度计算逻辑错误：代码中按单个坐标轴的距离平方计算加速度，完全不符合万有引力定律。万有引力加速度是基于三维空间中天体的直线距离，且是矢量，需要用单位矢量修正方向。
- 初始位置导致分母为0：天体2和天体3的x、y坐标完全相同，计算它们之间x、y轴的加速度时，分母为0，直接产生
inf，后续计算扩散为nan。 - 时间步长过大：100秒的步长对于天体运动来说太大，数值积分误差快速累积导致计算发散。
修复方案
1. 修正加速度计算逻辑
重新实现符合物理规律的加速度计算，先算三维距离，再通过单位矢量得到各轴加速度分量：
void calculate_acc() { // 计算天体1受到的加速度 long double dx12 = body_2_pos[0] - body_1_pos[0]; long double dy12 = body_2_pos[1] - body_1_pos[1]; long double dz12 = body_2_pos[2] - body_1_pos[2]; long double dist12_sq = dx12*dx12 + dy12*dy12 + dz12*dz12; long double dist12 = sqrt(dist12_sq); long double acc_mag12 = G * body_mass / dist12_sq; long double dx13 = body_3_pos[0] - body_1_pos[0]; long double dy13 = body_3_pos[1] - body_1_pos[1]; long double dz13 = body_3_pos[2] - body_1_pos[2]; long double dist13_sq = dx13*dx13 + dy13*dy13 + dz13*dz13; long double dist13 = sqrt(dist13_sq); long double acc_mag13 = G * body_mass / dist13_sq; acc_body_1[0] = acc_mag12 * dx12 / dist12 + acc_mag13 * dx13 / dist13; acc_body_1[1] = acc_mag12 * dy12 / dist12 + acc_mag13 * dy13 / dist13; acc_body_1[2] = acc_mag12 * dz12 / dist12 + acc_mag13 * dz13 / dist13; // 计算天体2受到的加速度 long double dx21 = body_1_pos[0] - body_2_pos[0]; long double dy21 = body_1_pos[1] - body_2_pos[1]; long double dz21 = body_1_pos[2] - body_2_pos[2]; long double dist21_sq = dx21*dx21 + dy21*dy21 + dz21*dz21; long double dist21 = sqrt(dist21_sq); long double acc_mag21 = G * body_mass / dist21_sq; long double dx23 = body_3_pos[0] - body_2_pos[0]; long double dy23 = body_3_pos[1] - body_2_pos[1]; long double dz23 = body_3_pos[2] - body_2_pos[2]; long double dist23_sq = dx23*dx23 + dy23*dy23 + dz23*dz23; long double dist23 = sqrt(dist23_sq); long double acc_mag23 = G * body_mass / dist23_sq; acc_body_2[0] = acc_mag21 * dx21 / dist21 + acc_mag23 * dx23 / dist23; acc_body_2[1] = acc_mag21 * dy21 / dist21 + acc_mag23 * dy23 / dist23; acc_body_2[2] = acc_mag21 * dz21 / dist21 + acc_mag23 * dz23 / dist23; // 计算天体3受到的加速度 long double dx31 = body_1_pos[0] - body_3_pos[0]; long double dy31 = body_1_pos[1] - body_3_pos[1]; long double dz31 = body_1_pos[2] - body_3_pos[2]; long double dist31_sq = dx31*dx31 + dy31*dy31 + dz31*dz31; long double dist31 = sqrt(dist31_sq); long double acc_mag31 = G * body_mass / dist31_sq; long double dx32 = body_2_pos[0] - body_3_pos[0]; long double dy32 = body_2_pos[1] - body_3_pos[1]; long double dz32 = body_2_pos[2] - body_3_pos[2]; long double dist32_sq = dx32*dx32 + dy32*dy32 + dz32*dz32; long double dist32 = sqrt(dist32_sq); long double acc_mag32 = G * body_mass / dist32_sq; acc_body_3[0] = acc_mag31 * dx31 / dist31 + acc_mag32 * dx32 / dist32; acc_body_3[1] = acc_mag31 * dy31 / dist31 + acc_mag32 * dy32 / dist32; acc_body_3[2] = acc_mag31 * dz31 / dist31 + acc_mag32 * dz32 / dist32; }
2. 修改初始位置避免分母为0
调整天体2和3的坐标，让它们的x/y轴存在差值：
long double array_2[3] ={10000000000, 10000000000, 10000000000}; long double array_3[3] = {10000000000 + 1000, 10000000000, -10000000000};
3. 缩小时间步长
将步长改为1秒，降低积分误差：
system.init_time_step(1);
4. 数值计算优化
- 用直接相乘替代
pow计算平方，精度更高且速度更快，比如dx*dx代替pow(dx,2)。 - 可以改用Verlet积分替代当前的欧拉积分，Verlet积分在长期模拟中稳定性更好，不易发散。
额外提示
指针解决不了这个问题，核心问题是物理公式错误和初始条件设置不当，和数据存储方式无关。
内容的提问来源于stack exchange，提问作者Ronny
相关产品推荐
相关产品推荐
