Java太阳系模拟器中月球无法绕地球公转的问题排查

问题背景 #

用Java开发基于真实数据的太阳系模拟器，核心逻辑基于牛顿万有引力公式：维护包含太阳的行星列表，每帧遍历列表计算每个行星受其他行星的合力，更新速度与位置。核心代码如下：

private Point2D getAttractionForce(Planet planet){
    double distanceX = planet.x-this.x;
    double distanceY = planet.y-this.y;
    double distance = Math.sqrt(Math.pow(distanceX, 2)+Math.pow(distanceY, 2));

    double force = SolarSystem.G*this.mass*planet.mass/Math.pow(distance, 2);
    double angle = Math.atan2(distanceY, distanceX);
    double forceX = force * Math.cos(angle);
    double forceY = force * Math.sin(angle);

    return new Point2D(forceX, forceY);
}

public void updatePosition(List<Planet> planets){
    double totalForceX = 0;
    double totalForceY = 0;

    for (Planet planet : planets){
        if (planet == this) continue;
            Point2D force = getAttractionForce(planet);
            totalForceX += force.getX();
            totalForceY += force.getY();
        }
    }

    this.xVelocity += totalForceX/this.mass*TIMESTEP;
    this.yVelocity += totalForceY/this.mass*TIMESTEP;

    this.x += this.xVelocity*TIMESTEP;
    this.y += this.yVelocity*TIMESTEP;
}

添加月球后，月球未按预期绕地球公转，反而被太阳捕获。调试数据显示：

Distance moon-sun: 1.504E11
Sun force for moon: 4.313523549126089E20
Distance moon-earth: 4.0E8
Earth force for moon: 1.831019119125E20

太阳对月球的引力是地球的4倍，调整月球初始速度后问题仍未解决。月球初始化代码：

Planet moon = new Planet(Color.GRAY, 4, 7.35e22, -(1.50e11+4e8), 0, 0, 1022);

问题根源与解决方案 #

1. 初始速度设置错误 #

月球的初始速度需叠加地球的公转速度：地球本身以约29783 m/s的速度绕太阳公转，月球的初始速度应为地球公转速度加上其相对地球的公转速度（约1022 m/s），而非单独设置为1022 m/s。

当前月球的初始yVelocity为1022，远低于地球的29783，导致月球相对于太阳的速度不足以维持与地球绑定的轨道，直接被太阳的强引力捕获。

修改后的月球初始化代码（匹配地球公转方向）：

// 地球的公转速度为29783 m/s，叠加月球相对地球的速度
Planet moon = new Planet(Color.GRAY, 4, 7.35e22, -(1.50e11+4e8), 0, 0, 29783 + 1022);

2. 欧拉积分误差累积 #

当前使用的欧拉积分方法误差较大，长期模拟会导致轨道偏移加剧，放大初始速度设置的问题。建议改用Verlet积分，该方法对轨道模拟的稳定性远高于欧拉积分：

1. 给Planet类添加previousX、previousY字段，初始化时与x、y赋值相同；
2. 修改updatePosition方法：

public void updatePosition(List<Planet> planets){
    double totalForceX = 0;
    double totalForceY = 0;

    for (Planet planet : planets){
        if (planet == this) continue;
        Point2D force = getAttractionForce(planet);
        totalForceX += force.getX();
        totalForceY += force.getY();
    }

    double ax = totalForceX / this.mass;
    double ay = totalForceY / this.mass;

    // Verlet积分：基于加速度更新位置，再推导速度
    double newX = 2 * this.x - this.previousX + ax * TIMESTEP * TIMESTEP;
    double newY = 2 * this.y - this.previousY + ay * TIMESTEP * TIMESTEP;

    this.xVelocity = (newX - this.previousX) / (2 * TIMESTEP);
    this.yVelocity = (newY - this.previousY) / (2 * TIMESTEP);

    // 更新位置记录
    this.previousX = this.x;
    this.previousY = this.y;
    this.x = newX;
    this.y = newY;
}

3. 时间步长优化 #

若TIMESTEP设置过大（如秒级以上），会进一步放大积分误差。建议缩小时间步长（如设为60秒），或采用自适应时间步长策略，平衡模拟精度与性能。

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