行星轨道的数值近似

Question

新的一天新的问题...

我按照埃里克说的方式做了，我创建了一个加速方法！

看起来像这样：

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 的方法接收一个带有位置的 Vector 2。现在我想使用这个向量，所以我创建了另一个方法，它应该 return 我的幅度作为标量（双精度）。为此，我创建了一个方法来计算向量的大小并将其平方。当我在我的 ForceAsVector-Method 中调用我的 MagnitudeOfForce-Method 时，我试图将一个负单位向量与一个标量相乘，并将其 return 作为 Vector2。

也许我在这里做得非常糟糕，但我正在努力理解 Eric 所做的一切，仍然需要你的一些帮助！

谢谢。

Answer 1

你犯了两个明显的错误和算法选择不当。

首先，您将月亮的位置建模为单个数字而不是 3-space 中的向量，因此您在这里建模的是简谐运动，例如 spring 受限于单一维度，而不是轨道。首先将位置和速度建模为 3-space 向量，而不是双数。

其次，你把引力的符号设为正，所以它是两个物体之间的排斥力，而不是 有吸引力 。力的符号必须在地球的方向。

第三，您对欧拉算法的实现似乎是正确的，但欧拉算法对于轨道问题的数值求解是一个糟糕的选择，因为它不保守；你很容易陷入月球获得或失去一点动量的情况，这会加起来，破坏你漂亮的椭圆轨道。

由于月球轨道是哈密顿量，所以改用辛算法；它旨在模拟保守系统。

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

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

https://en.wikipedia.org/wiki/Orbital_state_vectors

但是，对于简单的 2-body 系统，有更简单的方法。如果你想做的是像 Kerbal Space 程序这样的模拟，其中只有你正在运行的 body 影响你的轨道，而有多个卫星的行星是 "on rails"，你不会需要在每个时间单元上对系统进行仿真，计算出状态向量序列；您可以根据开普勒元素直接计算它们：

https://en.wikipedia.org/wiki/Orbital_elements

我们非常准确地知道月亮的六大元素；从那些你可以随时计算出月球在 3-space 中的位置，同样，假设没有任何东西干扰它的轨道。（实际上，月球的轨道是被太阳、地球不是球体、潮汐等等改变的。）

更新：

首先，因为这是课程作业你需要引用所有的来源，包括从互联网上获得帮助。你的导师必须知道什么工作是你的，什么工作是你让别人为你做的。

您问的是如何在二维空间中完成这项工作；这似乎是错误的，但不管怎样，按照课程作业说的去做。

我希望更多的初学者能学到的规则是：创建一个类型、方法或变量来解决你的问题。在这种情况下，我们希望表示复杂 value 的行为，因此它应该是 type，并且该类型应该是 值类型。在 C# 中，值类型是 struct。所以让我们这样做吧。

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})";

}

我们完成了向量。

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

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 来获得新的状态，如此重复。

作为初学者，仔细研究这些技巧。请注意我所做的：我创建了类型来表示概念，并创建了方法来表示这些概念上的操作。一旦你这样做了，程序就会变得非常清晰易读。

想想你将如何使用这些技术扩展这个程序；我们做了一个 hard-coded AccelerationOfTheMoon 方法，但是如果我们想计算多个物体的加速度怎么办？同样，创建一个类型来表示Body； body 的特征是 State 和质量。如果我们想解决 n-body 问题怎么办？我们的加速度计算必须将其他物体作为参数。等等。

行星轨道的数值近似

Numerical approximation of planetary orbit

c#

vector

newtons-method

numerical-methods

orbital-mechanics