Java LibGDX轨道模拟器速度等值异常 行星摆动无法正常公转

问题描述 #

正在使用LibGDX开发带初始速度的行星椭圆轨道模拟器，当前编写的天体运动计算逻辑如下：

public void move(float deltaTime, Planet planet) {
    float deltaX = planet.getPos().x - this.pos.x;
    float deltaY = planet.getPos().y - this.pos.y;
    float alpha = (float) Math.toDegrees(Math.atan2(deltaY, deltaX));

    float distance = (float) Math.sqrt(Math.pow(deltaX, 2) + Math.pow(deltaY, 2));
    float F = G * this.m * planet.getM() / distance*distance;
    this.force.x = F * MathUtils.cos(alpha);
    this.force.y = F * MathUtils.sin(alpha);

    this.vel.x += (this.force.x / this.m) * deltaTime;
    this.vel.y += (this.force.y / this.m) * deltaTime;

    this.pos.x += this.vel.x * deltaTime;
    this.pos.y += this.vel.y * deltaTime;
}

运行后行星出现无规律晃动，速度计算结果异常，完全无法形成稳定椭圆轨道。已核对万有引力物理原理无认知错误，判定为代码实现逻辑问题。
当前项目结构说明：

• 所有天体均为Planet类实例，暂未定义独立的恒星类
• Planet类仅包含Getter/Setter方法、与运动计算无关的render()方法，以及上述move()方法
  测试用两个天体的初始化参数如下：

planet1 = new Planet(30, 1, new Vector2(300, 300));
planet2 = new Planet(70, 332000, new Vector2(400, 400));

已知LibGDX圆形图形默认坐标原点为图形左下角而非圆心，已在构造函数中添加如下坐标偏移处理：

this.pos = pos;
this.pos.x -= r;
this.pos.y -= r;

代码错误点定位 #

共5处直接导致模拟失效的逻辑错误：

1. 万有引力公式运算优先级错误
   代码中float F = G * this.m * planet.getM() / distance*distance;的运算顺序为((G*m1*m2)/distance)*distance，最终结果为G*m1*m2，完全没有实现距离平方反比的计算，力的大小完全错误。正确写法应为G * this.m * planet.getM() / (distance * distance)。
2. 三角函数角度单位不匹配
   Math.atan2()返回值为弧度单位，代码中将其手动转为角度后，传入了要求入参为弧度的MathUtils.cos()/MathUtils.sin()方法，导致力的方向计算完全错误，这一步根本不需要做弧度转角度的处理，直接使用Math.atan2()返回的弧度值即可。
3. 物理坐标与渲染坐标混用
   物理计算需要使用天体质心坐标，但构造函数中直接将存储质心位置的pos变量减去了半径值，改成了渲染用的左下角坐标，导致两个天体的距离计算、力方向计算全部基于渲染点而非质心，结果必然出错。正确做法是pos始终存质心坐标，渲染圆形时再做坐标偏移。
4. 多体更新顺序错误
   当前逻辑是计算完一个天体的受力后立刻更新它的位置，再计算下一个天体的受力时就会用到已经更新过的、属于下一个时间步的位置，会导致系统能量不守恒，出现轨道晃动、发散的问题。正确做法是先遍历所有天体，计算出每个天体在当前时间步受到的总合力，暂存后再统一更新所有天体的速度、位置。
5. 初始参数缺失
   给出的初始化代码中没有给两个行星设置垂直于两天体连线的初始切向速度，初始速度为0的情况下两个天体只会在引力作用下直线相撞，根本不可能形成椭圆轨道。

修正后参考实现 #

构造函数修正（分离物理坐标和渲染坐标） #

// pos参数传入质心初始位置，直接存储为物理计算用坐标
this.pos = pos;
this.radius = r;
this.mass = m;
this.vel = new Vector2(0,0);
this.force = new Vector2(0,0);

运动逻辑修正（先算力再统一更新） #

// 每帧先重置所有天体的受力为0
public void resetForce() {
    this.force.set(0, 0);
}

// 计算当前天体受到来自otherPlanet的引力，累加受力，不更新位置
public void addForceFrom(Planet otherPlanet) {
    float deltaX = otherPlanet.pos.x - this.pos.x;
    float deltaY = otherPlanet.pos.y - this.pos.y;
    float distanceSq = deltaX*deltaX + deltaY*deltaY;
    float distance = (float) Math.sqrt(distanceSq);
    // 避免距离过近除0，加最小距离阈值
    if (distance < 1f) distance = 1f;
    float F = G * this.mass * otherPlanet.mass / distanceSq;
    // 直接用向量归一化计算力的分量，不需要转角度，避免三角函数单位错误
    this.force.x += F * deltaX / distance;
    this.force.y += F * deltaY / distance;
}

// 所有天体受力计算完成后，统一调用更新速度和位置
public void update(float deltaTime) {
    this.vel.x += (this.force.x / this.mass) * deltaTime;
    this.vel.y += (this.force.y / this.mass) * deltaTime;
    this.pos.x += this.vel.x * deltaTime;
    this.pos.y += this.vel.y * deltaTime;
}

// 渲染时再做坐标偏移，把质心坐标转成圆形左下角坐标
public void render(ShapeRenderer shapeRenderer) {
    shapeRenderer.circle(this.pos.x - radius, this.pos.y - radius, radius);
}

主循环调用逻辑 #

// 每帧执行顺序
// 1. 重置所有行星受力
for (Planet p : planets) {
    p.resetForce();
}
// 2. 两两计算引力累加，符合牛顿第三定律
for (int i = 0; i < planets.size; i++) {
    for (int j = i+1; j < planets.size; j++) {
        Planet p1 = planets.get(i);
        Planet p2 = planets.get(j);
        p1.addForceFrom(p2);
        p2.addForceFrom(p1);
    }
}
// 3. 统一更新所有行星位置速度
for (Planet p : planets) {
    p.update(deltaTime);
}
// 4. 渲染
for (Planet p : planets) {
    p.render(shapeRenderer);
}

初始速度设置参考 #

要形成稳定轨道，需要给小质量天体设置切向初速度，例如两个天体初始距离为√2*100≈141像素，可给质量1的小行星设置垂直于连线方向的初速度，大小根据实际G值调试即可：

// 大质量行星放在(400,400)，初始速度设为0（因为质量远大于小行星，近似不动）
planet2 = new Planet(70, 332000, new Vector2(400, 400));
// 小行星放在(300,300)，连线方向是左上到右下，切向速度设为(3, -3)方向
planet1 = new Planet(30, 1, new Vector2(300, 300));
planet1.vel.set(3f, -3f);

注意：如果轨道不稳定出现缓慢外扩/内旋，可适当调小deltaTime，上述实现使用的半隐式欧拉法已经能满足基础轨道模拟的精度要求，若需要更高精度可换用RK4积分方法。

内容的提问来源于stack exchange，提问作者eisa.exe