圆形轨道中重力的矢量方向

Question

我目前在使用C#进行一个项目，其中探索行星引力，这是一个很难理解的主题，但我喜欢挑战。我一直在研究牛顿定律和开普勒定律，但有一件事我无法解决，那就是如何得出正确的引力方向。

在我的示例中，我只有两个物体。一个是卫星，一个是行星。这是为了简化它，以便我能够理解 - 但我的计划是有多个物体动态地相互影响，并希望最终得到一个比较真实的多体系统。

当你有一个轨道时，那么卫星会受到引力作用，当然是朝向行星的方向，但该方向并不是恒定的。为了更好地解释我的问题，我将尝试使用一个例子：

假设我们有一个卫星，速度为50米/秒，并以10米/秒²的速度加速向行星移动，半径为100米。（所有数字仅为理论值）如果我们假设帧速率为1，则经过一秒钟后，物体将向前移动50个单位并向下移动10个单位。

由于卫星在一个帧中移动了多个单位，大约为半径的50％，因此重力方向在这个帧中已经发生了很大的变化，但施加的力只是“向下”。这会导致非常大的误差，特别是如果物体移动了半径的一大部分。

在我们的示例中，我们可能需要将引力方向基于当前位置和本帧结束时的位置之间的平均值。

如何计算这个问题呢？

我对三角学有基础的理解，但主要集中在三角形上。请假设我很蠢，因为与你们任何人相比，我可能都是。

（我之前提出过一个问题，但最终删除了它，因为它引起了一些敌意，而且基本上没有用语言表述清楚，并且太过笼统 - 它并不是一个具体的问题。我希望这次更好一些。如果还不行，请告诉我，我来这里是为了学习 :)）

仅供参考，这是我目前用于运动的函数：

foreach (ExtTerBody OtherObject in UniverseController.CurrentUniverse.ExterTerBodies.Where(x => x != this))
{
    double massOther = OtherObject.Mass;

    double R = Vector2Math.Distance(Position, OtherObject.Position);

    double V = (massOther) / Math.Pow(R,2) * UniverseController.DeltaTime;

    Vector2 NonNormTwo = (OtherObject.Position - Position).Normalized() * V;

    Vector2 NonNormDir = Velocity + NonNormTwo;
    Velocity = NonNormDir;

    Position += Velocity * Time.DeltaTime;
}

如果我表达得不好，请让我重新表达 - 英语不是我的母语，特定主题可能很难表达，当您不知道正确的技术术语时。 :)

我有一种预感，这在开普勒第二定律中有所涵盖，但如果是这样，那么我不确定如何使用它，因为我不完全理解他的定律。

谢谢您抽出时间 - 这意味着很多！

（还有，如果有人看到我功能上的多个错误，请指出！）

- Casper Bang

---

这是所有离散时间步积分器的一个问题，不仅仅局限于模拟重力。你可能需要在这里稍微了解一下：http://gafferongames.com/game-physics/integration-basics/ :) - Magnus Hoff

@MagnusHoff 这真是太棒了。谢谢 - 我会阅读并在尝试使用它来解决我的问题后回报！ :) - Casper Bang

我再怎么强调也不为过，这个问题非常重要：天文模拟中正确的积分方法是什么？ - user2005819

Answer 1

我目前在一个C#项目中，通过模拟行星引力来进行实验。这是一种同时学习模拟技术、编程和物理的有趣方式。

有一件事我无法解决，那就是如何获得正确的引力方向。我假设你不想模拟相对论引力。地球不是绕着太阳转，而是绕着太阳八分钟前所在的位置转。纠正引力不是瞬时的这一事实可能很困难。（更新：根据评论，这是不正确的。我只学了牛顿动力学的第二年，对张量微积分的理解非常模糊。）

在早期阶段，最好假设引力是瞬时的，行星是质点，所有质量都集中在中心。引力作用的矢量是从一个点到另一个点的直线。

假设我们有一个卫星以50m/s的速度运动……如果我们认为帧率是每秒一帧，那么经过一秒钟，物体将向右移动50个单位并向下移动10个单位。

让我们更清楚地说明一下。力等于质量乘以加速度。你计算出两个物体之间的力。你知道它们的质量，因此你现在知道每个物体的加速度。每个物体都有一个位置和一个速度。加速度改变速度，速度改变位置。因此，如果粒子开始时具有50m/s向左和0m/s向下的速度，然后你施加一个加速度为10m/s²向下的力，那么我们可以计算出速度的变化，然后是位置的变化。正如你所指出的，在那一秒钟结束时，位置和速度的变化都将与它们现有的大小相比发生巨大变化。

随着卫星在一个帧内移动多个单位和约50%的半径，引力方向在这个帧期间已经发生了很大的偏移，但施加的力只是“向下”的。这会产生很大的误差，特别是如果物体移动了大比例的半径。

正确。问题在于帧率太低，无法正确模拟你描述的相互作用。如果物体的方向变化如此迅速，你需要以十分之一、百分之一或千分之一秒为单位运行模拟。时间步长的大小通常被称为模拟的“delta t”，而你的模拟步长太大了。

对于行星体来说，你现在正在试图通过每几个月模拟其位置并假设其在此期间内沿直线运动来模拟地球。你需要每隔几分钟模拟一次其位置，而不是每隔几个月。

在我们的例子中，我们可能需要将引力方向基于当前位置和本帧末位置之间的平均值。

你可以这样做，但更容易的方法是减小计算中的“delta t”值。这样，在帧的开始和结束时方向之间的差异就会更小。

解决了这个问题后，你可以使用更多的技巧。例如，当位置在帧之间发生太大变化时，可以“检测”并返回以较小时间步长重新计算。如果位置几乎没有变化，则增加时间步长。

一旦解决了这个问题，就可以在物理模拟中使用很多更高级的技术，但我建议首先确立基本的时间步进。更高级的技术本质上是对你的想法“智能插值改变时间步长”的变体，你在这方面已经走得不错了，但先学会基础再进阶比较好。

- Eric Lippert

---

谢谢你的回答。正如你所看到的，“UniverseController.DeltaTime;”已经是我的计算的一部分，它是上一帧所需的时间乘以一个时间比例。我知道这是个问题，但我想能够以更高的速度移动我的卫星/粒子而不会失去太多精度。必须说，我并不确定这是唯一的问题，可能还有计算和实现方面的问题，但这似乎是错误的一个来源，我可以进行修正。当然，如果其他问题也能得到解决，那就更好了。 - Casper Bang

@CDK：假设模拟已经按比例缩放，使得万有引力常数为1，则计算似乎是正确的。加速度为m1 x m2 / r ^ 2，因此具有质量m1的粒子在时间跨度delta t内的delta v是dt * m1 x m2 / r ^ 2 / m1，m1相互抵消，您就得到了正确的delta v计算。只需将delta t变得更小，您的模拟就会更准确。 - Eric Lippert

@CDK 你可以用RK4得到一个相当不错的近似值，但是对于任何合理的时间步长来说，欧拉方法并不可靠（它不够准确）。如果你能想出一个一阶微分方程状态向量，并将其应用于RK4，那么这个问题就有了相当优雅的解决方案，但我不确定你的数学背景有多好。 - sapi

“地球不是绕着太阳轨道运行，而是绕着太阳八分钟前所在的位置运行。”这种说法非常错误，早在1805年就被拉普拉斯发现了。爱因斯坦110年后才知道这一点；他的广义相对论并没有这样说。完整的广义相对论公式会将其转化为偏微分方程问题，而一阶参数化后牛顿近似则仍处于常微分方程领域，并且对于太阳系工作得非常好；请参见《太阳、月亮和行星的轨道历书》第3页的8-1式方程。 - David Hammen

只是提醒一下，rk4在长时间内仍然会不稳定，因为它不能保持引力系统的能量。您需要一个辛积分器。第一阶方法仍然需要您的积分器的一半是隐式的，但它将永远保持稳定，不像标准的欧拉或rk4。 - Godric Seer

Answer 2

我将从一种技术开始介绍，它几乎和你一直在使用的欧拉-克罗默积分一样简单，但精度更高。这就是跳跃法（Leapfrog Technique）。其思想非常简单：位置和速度彼此保持半个时间步长。

初始状态下，位置和速度位于时间t0。为了获得半个时间步长的偏移量，您需要一个特殊情况来处理第一步，其中速度使用间隔开始时的加速度提前半个时间步长，然后位置向前移动一个完整步骤。在这个第一次特殊情况之后，代码就像您的欧拉-克罗默积分器一样工作。

伪代码如下：

void calculate_accel (orbiting_body_collection, central_body) {
    foreach (orbiting_body : orbiting_body_collection) {
        delta_pos = central_body.pos - orbiting_body.pos;
        orbiting_body.acc =
            (central_body.mu / pow(delta_pos.magnitude(),3)) * delta_pos;
    }
}

void leapfrog_step (orbiting_body_collection, central_body, delta_t) {
    static bool initialized = false;
    calculate_accel (orbiting_body_collection, central_body);
    if (! initialized) {
        initialized = true;
        foreach orbiting_body {
            orbiting_body.vel += orbiting_body.acc*delta_t/2.0;
            orbiting_body.pos += orbiting_body.vel*delta_t;
        }
    }
    else {
        foreach orbiting_body {
            orbiting_body.vel += orbiting_body.acc*delta_t;
            orbiting_body.pos += orbiting_body.vel*delta_t;
        }
    }
}

请注意，我已将加速度作为每个轨道体的字段添加。这是一个临时的步骤，目的是使算法与你的算法类似。请注意，我将加速度的计算移到了它自己的独立函数中。这不是一个临时的步骤。这是进一步提高集成技术所必需的第一步。

下一步关键的步骤是撤销对加速度的临时添加。加速度应该属于积分器，而不是物体本身。另一方面，加速度的计算属于问题空间，而不是积分器。您可能想添加相对论修正、太阳辐射压力或行星之间的重力相互作用。积分器不应该知道这些加速度是如何计算出来的。函数“calculate_accels”是一个黑盒子，由积分器调用。

不同的积分器对计算加速度的时间概念非常不同。有些会存储最近加速度的历史记录，有些需要额外的工作空间来计算某种平均加速度。有些也会使用速度(保留历史记录，拥有一些速度工作空间)。一些更先进的集成技术在内部使用许多技术，从一个技术切换到另一个技术，以提供最佳的精度和CPU使用率平衡。如果您想模拟太阳系，您需要一个非常精确的积分器。(而且你需要远离浮点数。即使是双倍精度对于高精度的太阳系积分也不够好。使用浮点数，没有必要超越RK4，甚至可能不需要leapfrog。)

正确地分离属于谁(积分器与问题空间)使得可以精细化问题领域(添加相对论等)，并且易于切换集成技术，以便您可以评估一种技术与另一种技术。

- David Hammen

---

这可能正是我需要的，我会记下来。但目前它对我来说太高级了，无法理解。顺便说一句，在计算实际运动时，我已经只使用双精度浮点数了，但什么更好呢？我会尝试多读几遍你的东西，但现在，我提供的答案还可以，它消除了任何大时间步进的可能性，因此相当准确。但我会朝着你的想法努力，并看看它是否能胜任工作。我的计划是添加火箭和更多功能，直到这个项目让我感到无聊为止。谢谢你的时间！ :) - Casper Bang

从您当前使用的Euler-Cromer技术切换到Leapfrog相对容易，并且增加准确性的回报非常巨大。您可以通过仅进行微小的更改在当前架构中进行此切换。在这里显示的代码使用浮点数而不是双精度。当涉及建模轨道力学时，浮点数很快就会崩溃。从浮点数切换到双精度可获得很大的回报，而几乎不会影响性能。（续） - David Hammen

关于更高的精度：许多机器不支持本地化。在C/C++中的long double可能只有80位宽（例如，英特尔盒子）。为了避免即使使用最好的积分器也会完全失去精度，模拟几千年甚至更长时间的太阳系需要至少四倍精度（128位浮点数）。这意味着在不支持四倍精度的计算机上进行软件仿真，这反过来又意味着极大的性能损失。在您的应用程序中没有理由去那里。双精度应该对您而言足够好。 - David Hammen

我真的不确定如何实现上述内容。我可以在哪里联系您以了解有关此概念的问题？无论如何，如果不行： 我上面展示的代码使用双精度浮点数.. "double V = (massOther) / Math.Pow(R,2) * Time.DeltaTime;" 等等... 为什么第一帧要使用一半的加速度？"orbiting_body.vel += orbiting_body.acc*delta_t/2.0;" 我得到的和你提出的主要区别是什么？我看不出有什么大的区别.. 抱歉。 我会尝试在需要的地方添加小数.. 或者我应该使用其他数据类型吗？谢谢迄今为止的帮助.. - Casper Bang

Answer 3

所以我找到了一个解决方案，可能不是最聪明的，但它有效，并且在阅读了Eric的答案和Marcus的评论后，这个想法很快就浮现出来了，你可以说它是两者的结合：

这是新代码：

foreach (ExtTerBody OtherObject in UniverseController.CurrentUniverse.ExterTerBodies.Where(x =>  x != this))
{
double massOther = OtherObject.Mass;

double R = Vector2Math.Distance(Position, OtherObject.Position);

double V = (massOther) / Math.Pow(R,2) * Time.DeltaTime;

float VRmod = (float)Math.Round(V/(R*0.001), 0, MidpointRounding.AwayFromZero);
if(V > R*0.01f)
{
    for (int x = 0; x < VRmod; x++)
    {
        EulerMovement(OtherObject, Time.DeltaTime / VRmod);
    }
}
else
    EulerMovement(OtherObject, Time.DeltaTime);

}

public void EulerMovement(ExtTerBody OtherObject, float deltaTime)
    {

            double massOther = OtherObject.Mass;

            double R = Vector2Math.Distance(Position, OtherObject.Position);

            double V = (massOther) / Math.Pow(R, 2) * deltaTime;

            Vector2 NonNormTwo = (OtherObject.Position - Position).Normalized() * V;

            Vector2 NonNormDir = Velocity + NonNormTwo;
            Velocity = NonNormDir;



            //Debug.WriteLine("Velocity=" + Velocity);
            Position += Velocity * deltaTime;
    }

解释一下：

我得出的结论是，如果卫星在一个框架中的速度过快，那么为什么不将其分成多个框架呢？这就是现在的“它”所做的。

当卫星的速度超过当前半径的1％时，它会将计算分成多个部分，使其更精确。当处理高速度时，帧率当然会降低，但对于这样的项目来说，这是可以接受的。

不同的解决方案仍然非常受欢迎。我可能会调整触发量，但最重要的是它能够正常工作，然后我可以担心如何使其更加平滑！

感谢所有看过并帮助我自己得出结论的人们！ :) 人们能够像这样提供帮助真是太棒了！

- Casper Bang

---

请不要将更新发布为新答案；只需以答案的形式发布答案。编辑您的原始帖子以添加新信息。 - Eric Lippert

实际上，回答问题的更新应该作为答案发布，而不是作为原始帖子的编辑。CDK做得很正确。 - Brian