行星轨道的数值逼近

Question

新的一天新的问题...

我按照埃里克（Eric）所说的做，并为Acceleration创建了一种方法！

看起来像这样：

static Vector2 AccelerationOfTheMoon(Vector2 position)
    {
        double Mm = 7.349 * Math.Pow(10, 22);

        double MagnitudeOfForce()
        {
            double MagF = position.LengthSquared();
            return MagF;
        }

        Vector2 ForceAsVector()
        {
            return position.NegativeUnitVector() * MagnitudeOfForce();
        }

        Vector2 acceleration = ForceAsVector() / Mm;


    }

以我的方式看，我的名为AccelerationOfTheMoon的方法会收到带有位置的向量2。现在，我想使用此向量，因此我创建了另一个方法，该方法应将我的Magnitude作为标量（双精度）返回。为此，我创建了一种方法来计算向量的幅度并将其平方。当我在ForceAsVector-Method中调用MagnitudeOfForce-Method时，我正在尝试将负单位矢量与标量相乘，并将其作为Vector2返回。

也许我在这里做得很糟糕，但是我正在努力了解Eric所做的一切，仍然需要您的帮助！

谢谢。

Answer 1

您犯了两个明显的错误，并且选择了错误的算法。

首先，您已经将月球的位置建模为单个数字，而不是在3空间中的向量，所以这里要建模的是简单的谐波运动，例如约束在一个维度上的弹簧，而不是轨道。首先将位置和速度建模为3空间向量，而不是两倍。

其次，您已经使引力的符号为正，因此它是两个物体之间的排斥力，而不是吸引力。力的符号必须指向地球的方向。

第三，您对Euler算法的实现似乎是正确的，但是Euler对于数值求解轨道问题不是一个好的选择，因为它并不保守。您可以轻松进入月球获得或失去一点动量的情况，这些情况加起来并破坏了您的椭圆形轨道。

由于月球的轨道是哈密顿量，因此请使用辛算法。它旨在模拟保守系统。

https://zh.wikipedia.org/wiki/Symplectic_integrator

此方法和您的欧拉方法从根本上讲是关于寻找状态向量的：

https://zh.wikipedia.org/wiki/轨道状态_矢量

但是，对于简单的2体系统，有更简便的方法。如果您要做的是像Kerbal Space Program这样的模拟，其中只有您所运行的物体会影响您的轨道，而有多个卫星的行星都在“轨道上”，那么您无需在每个时间单位上都对系统进行模拟计算状态向量的序列；您可以根据开普勒元素直接计算它们：

https://zh.wikipedia.org/wiki/轨道元素

我们高度了解月球的六个元素；再次假设那些都没有扰动其轨道的情况，您可以随时从中计算出其在月球3空间中的位置。（实际上，月球的轨道因太阳，地球不是球体，由于潮汐等而改变。）

更新：

首先，由于这是课程工作，因此您必须引用所有资源，包括从互联网获得帮助。您的讲师必须知道您的工作是什么，以及别人为您所做的工作。

您问过如何在两个方面进行这项工作；这似乎是错误的，但是无论如何，要按照课程的工作去做。

我希望能教给更多的初学者的规则是：创建一种类型，方法或变量来解决您的问题。在这种情况下，我们希望表示一个复杂值的行为，因此它应该是一个类型，而该类型应该是一个值类型。值类型struct在C＃中。因此，让我们这样做。

struct Vector2
{
    double X { get; }
    double Y { get; }
    public Vector2(double x, double y)
    {
        this.X = x;
        this.Y = y;
    }

请注意，向量是不变的，就像数字一样。您永远不会变异向量。当您需要一个新矢量时，就可以制作一个新矢量。

我们需要对向量执行哪些操作？向量加法，标量乘法和标量除法只是花式乘法。让我们实现这些：

    public static Vector2 operator +(Vector2 a, Vector2 b) => 
      new Vector2(a.X + b.X, a.Y + b.Y);
    public static Vector2 operator -(Vector2 a, Vector2 b) => 
      new Vector2(a.X - b.X, a.Y - b.Y);
    public static Vector2 operator *(Vector2 a, double b) => 
      new Vector2(a.X * b, a.Y * b);
    public static Vector2 operator /(Vector2 a, double b) => 
      a * (1.0 / b);

我们也可以按其他顺序进行乘法，因此让我们实现它：

    public static Vector2 operator *(double b, Vector2 a) => 
      a * b;

将向量设为负数等于将其乘以-1：

    public static Vector2 operator -(Vector2 a) => a * -1.0;

并帮助我们调试：

    public override string ToString() => $"({this.X},{this.Y})";

}

向量已经完成。

我们正在尝试模拟轨道状态参数的演化，因此再次创建一个type。状态参数是什么？位置和速度：

struct State
{
    Vector2 Position { get; }
    Vector2 Velocity { get; }
    public State(Vector2 position, Vector2 velocity)
    {
        this.Position = position;
        this.Velocity = velocity;
    }

现在我们来了解核心算法。我们有什么？状态和加速度。我们想要什么？新状态。制作方法：

    public State Euler(Vector2 acceleration, double step)
    {
        Vector2 newVelocity = this.Velocity + acceleration * step;
        Vector2 newPosition = this.Position + this.Velocity * step;
        return new State(newPosition, newVelocity);
    }
}

超。还剩什么？ 我们需要解决加速问题。制作方法：

static Vector2 AccelerationOfTheMoon(Vector2 position) 
{
  // You do this. Acceleration is force divided by mass,
  // force is a vector, mass is a scalar. What is the force
  // given the position? DO NOT MAKE A SIGN MISTAKE AGAIN.
}

现在，您拥有了所需的所有部件。从初始状态开始，您可以重复调用AccelerationOfTheMoon获取新的加速度，然后调用Euler获取新的状态并重复。

作为初学者，请仔细研究这些技术。注意我做了什么：我创建了代表概念的类型，并创建了代表这些概念上的操作的方法。完成此操作后，该程序将变得非常清晰易读。

考虑使用这些技术如何扩展该程序；我们做了一个硬编码的AccelerationOfTheMoon方法，但是如果我们想计算几个物体的加速度怎么办？再次，使类型表示Body；身体的特征是State质量。如果我们想解决n体问题怎么办？我们的加速度计算必须将其他物体作为参数。等等。