在 C# 中计算自定义太阳系

Question

我正在尝试为一个研究项目创建一个定制的太阳系，以测试不同的数值方法并了解它们的差异。

我们正在建造自己的太阳系，而不是尝试模拟地球/太阳，因为这看起来更困难。

现在我们有太阳 M1 和地球 M2。

• M1 = 3500
• M2=500
• M1 到 M2 的距离 = 200
• 我们使用不同的 G = 0.2

当我们计算重力时，F = 8.75 或 -8.75 (F = GM1M2/r^2)

当我们计算地球绕地球公转的恒定速度时，我们发现 v = 2 (v = sqrt(G(M1+M2)/r))

为了计算力的 vec2，我们使用以下代码

    static double GravitationalConstant = 0.2f;//6.674e-11;
    static public Vector2f GravitationalForce(SolarObject onObj, SolarObject fromObj)
    {
        Vector2f result = fromObj.OldPosition - onObj.OldPosition;

        float distance = result.Lenght();
        Console.WriteLine(distance);
        result = new Vector2f(result.X / distance, result.Y / distance); // Normalized, but I've the lenght al ready so

        result *= (float)((GravitationalConstant * onObj.Mass * fromObj.Mass) / (distance * distance));

        return result;
    }

为了更新位置/速度，我们使用这个

public void Update(float dt, List<SolarObject> objects)
    {
        foreach(SolarObject s in objects.Skip(1))
        {
            Vector2f f = Utility.GravitationalForce(s, objects[0]);

            Vector2f a = f / s.Mass;

            s.OldPosition  = s.NewPosition;
            s.NewPosition += dt * s.Velocity
            s.Velocity += dt * a;
        }
    }

该物体确实绕太阳飞行，但它根本不是轨道。M1/M2 之间的距离不是恒定的 <-> 力也不总是等于 8.75f。我们知道欧拉有一个“误差”，但这似乎很大，因为即使没有绕一圈，距离也已经增加了 25%。所以一定是某个地方有 bug。

Answer 1

不幸的是，这种行为是欧拉积分所固有的——你总是会在某种程度上超出弯曲路径。这种影响可以通过使用更小的时间步来抑制，即使没有doubles 也能起作用：

• dt = 0.01：

  

• dt = 0.001：

  

• dt = 0.0001：

  

• dt = 0.00001：

  

正如您所看到的，欧拉方法的准确性随着时间步长的减小而提高。内行星（较小的轨道半径=较大的曲率=更多的超调）开始更一致地遵循其投影轨道（绿色）。超调仅在最内层行星 处可见dt = 0.0001，而在 处则完全不可见dt = 0.00001。

为了改进欧拉方法而不必求助于极其小的时间步长，可以使用例如龙格-库塔积分（四阶变体很流行）。

此外，轨道速度应该是v = sqrt(G*Msun/r))而不是v = sqrt(G(M1+M2)/r))，尽管对于大恒星的限制，这不会造成太大的问题。

（如果您想要我的测试代码，请告诉我 - 尽管它写得非常糟糕，并且核心功能与您的相同。）