有哪些算法可以让我模拟行星物理学？

Question

有哪些算法可以让我模拟行星物理学？

algorithmsimulationphysicsastronomy

28

我对制作一个“太阳系”模拟器很感兴趣，可以模拟行星和恒星的旋转和引力力量。

我希望能够模拟我们的太阳系，并在不同的速度下进行模拟（例如，观察地球和其他行星在几天、几年等时间内绕太阳旋转）。我想能够添加行星并更改它们的质量等参数，以查看它们对系统的影响。

有没有资源可以指导我编写这种模拟器？是否有任何现有的物理引擎专门设计用于此目的？

- Gaius

6

万有引力定律公式为：F = G.m1.m2.r^2，其中G为引力常数，m1和m2为两个物体的质量，r为它们之间的距离。 - skaffman

14

@skaffman：r^(-2) ;) - Stefano Borini

@skaffman 我想要这种公式，但是不包括碰撞，用于模拟“流体”。我对数学/物理并不是很了解，所以如果你能给我一些简单的帮助，我会非常感激。 - Llamageddon

能否在大小和质量方面制作逼真的n体太阳系模拟？请参见此处。 - Spektre

12个回答

8

您需要了解和理解牛顿万有引力定律和开普勒行星运动定律。这两个定律都很简单，我相信你在高中时应该听说过它们，如果没有学习过，也应该有所耳闻。最后，如果您想让模拟器尽可能准确，您应该熟悉n-体问题。

您应该从简单的开始。尝试制作一个绕着太阳旋转的太阳对象和一个地球对象。那将为您提供非常坚实的起点，并且从那里扩展相对容易。一个行星对象看起来应该是这样的：

Class Planet {
  float x;
  float y;
  float z; // If you want to work in 3D
  double velocity;
  int mass;
}

只要记住 F = MA，其余的只是无聊的数学:P

- David Titarenco

3

你可能想考虑使用极坐标。在轨道力学中，它们通常更容易使用。 - MSalters

1

没错，但你需要在屏幕上绘制行星（我猜），所以最好在笛卡尔平面上进行初始计算。 - David Titarenco

在直角坐标系和极坐标系之间进行转换是微不足道的，因此可以根据需要方便地进行计算并进行必要的转换。 - David Thornley

6

这是一篇关于N体问题的优秀教程。

http://www.artcompsci.org/#msa

它使用Ruby编写，但很容易映射到其他语言中。其中包括了一些常见的积分方法：Forward-Euler、Leapfrog和Hermite。

- Ade Miller

4

您可能想要看一下Celestia，这是一个免费的太空模拟器。我相信您可以使用它来创建虚构的太阳系，并且它是开源的。

- High Performance Mark

2

你只需要实现正确的微分方程（开普勒定律）并使用龙格-库塔方法即可。（至少对我有用，但可能有更好的方法）。有许多这样的模拟器在线上。以下是一个简单的实现了500行C代码的模拟器（运动算法要少得多）：http://astro.berkeley.edu/~dperley/programs/ssms.html。还可以查看以下内容：
http://en.wikipedia.org/wiki/Kepler_problem
http://en.wikipedia.org/wiki/Two-body_problem
http://en.wikipedia.org/wiki/N-body_problem。

- Luka Rahne

7

不要使用龙格-库塔(RK4)，而是使用速度Verlet(或更高阶的辛方法，如Lobatto IIIA-IIIB)。RK4虽然高阶，但缺乏后者保存结构的特性，使用它会导致你的行星最终漂出太空或撞向太阳。 - user168715

你也应该发表一个答案，未知先生！ - Gaius

1

在高中时，我们研究了几种数值求导方法。Leapfrog（速度Verlet）的特点是它不能保持能量不变。（微分方程实际上是开普勒定律）。从我所检查的情况来看，龙格-库塔方法可以保持能量不变（其他方法有所谓的能量漂移），我实现了大约5种方法，如果有人想看我的第一次编程尝试，应该还有一些Matlab文件在我的备份中。 - Luka Rahne

2

这似乎与维基百科（http://en.wikipedia.org/wiki/Leapfrog_integration）不符：Leapfrog积分的第二个优点是其辛性质，这意味着它保留动力系统的（稍作修改的）能量。当计算轨道动力学时，这特别有用，因为其他积分方案（如Runge-Kutta方法）不保留能量并允许系统随时间大幅漂移。你有参考资料吗？ - idontgetoutmuch

2

在物理学中，这被称为N体问题。它很有名，因为你无法手动解决超过三个行星的系统。幸运的是，你可以用计算机轻松获得近似解决方案。

一个关于从头开始编写此代码的好文章可以在这里找到。

然而，我觉得需要警告一句话。你可能得不到你期望的结果。如果你想看看：

行星的质量如何影响其围绕太阳的轨道速度，那很酷。你会看到这一点。
不同的行星相互作用，你会失望的。

问题就在这里。 是的，现代天文学家关心土星的质量如何改变地球围绕太阳的轨道。但这只是一个非常小的影响。如果你要绘制行星围绕太阳的路径，其他太阳系中的行星几乎不重要。太阳太大了，它将淹没所有其他引力。唯一的例外是：

如果你的行星轨道非常椭圆，这将导致行星之间可能更加接近，因此它们会相互作用。

如果你的行星与太阳的距离几乎相同，它们之间也会相互作用。

如果你让你的行星变得非常巨大，以至于它们在外太阳系中与太阳竞争引力，那么它们也会相互作用。

需要明确的是，如果你创建一个真实的太阳系，是可以计算出一些行星之间的相互作用的，但是这些相互作用对于肉眼来说并不显著。

不过，你可以尝试一下，看看会怎样！

- john_science

1

用于模拟行星物理学的算法。

这是我在Android应用程序中实现Keppler部分的方法。主要部分可以在我的网站上下载整个源代码：http://www.barrythomas.co.uk/keppler.html

这是我绘制行星在轨道上“下一个”位置的方法。想象一下，像在一个圆上每次向前走一度那样，按照与你试图跟踪的行星相同周期的圆形步进。在这个方法之外，我使用一个全局double作为步数计数器——称为dTime，其中包含旋转的度数。

传递给该方法的关键参数是dEccentricty，dScalar（一个缩放因子，使轨道全部适合于显示），dYear（地球年份内的轨道持续时间）以及将轨道定位到近日点位于表盘正确位置的dLongPeri——近日点的经度。

drawPlanet：

public void drawPlanet (double dEccentricity, double dScalar, double dYear, Canvas canvas, Paint paint, 
            String sName, Bitmap bmp, double dLongPeri)
{
        double dE, dr, dv, dSatX, dSatY, dSatXCorrected, dSatYCorrected;
        float fX, fY;
        int iSunXOffset = getWidth() / 2;
        int iSunYOffset = getHeight() / 2;

        // get the value of E from the angle travelled in this 'tick'

        dE = getE (dTime * (1 / dYear), dEccentricity);

        // get r: the length of 'radius' vector

        dr = getRfromE (dE, dEccentricity, dScalar);

        // calculate v - the true anomaly

        dv = 2 * Math.atan (
                Math.sqrt((1 + dEccentricity) / (1 - dEccentricity))
                *
                Math.tan(dE / 2)
                ); 

        // get X and Y coords based on the origin

        dSatX = dr / Math.sin(Math.PI / 2) * Math.sin(dv);
        dSatY = Math.sin((Math.PI / 2) - dv) * (dSatX / Math.sin(dv));

        // now correct for Longitude of Perihelion for this planet

        dSatXCorrected = dSatX * (float)Math.cos (Math.toRadians(dLongPeri)) - 
            dSatY * (float)Math.sin(Math.toRadians(dLongPeri));
        dSatYCorrected = dSatX * (float)Math.sin (Math.toRadians(dLongPeri)) + 
            dSatY * (float)Math.cos(Math.toRadians(dLongPeri));

        // offset the origin to nearer the centre of the display

        fX = (float)dSatXCorrected + (float)iSunXOffset;
        fY = (float)dSatYCorrected + (float)iSunYOffset;

        if (bDrawOrbits)
            {
            // draw the path of the orbit travelled
            paint.setColor(Color.WHITE);
            paint.setStyle(Paint.Style.STROKE);
            paint.setAntiAlias(true);

            // get the size of the rect which encloses the elliptical orbit

            dE = getE (0.0, dEccentricity);
            dr = getRfromE (dE, dEccentricity, dScalar);
            rectOval.bottom = (float)dr;
            dE = getE (180.0, dEccentricity);
            dr = getRfromE (dE, dEccentricity, dScalar);
            rectOval.top = (float)(0 - dr);

            // calculate minor axis from major axis and eccentricity
            // http://www.1728.org/ellipse.htm

            double dMajor = rectOval.bottom - rectOval.top;
            double dMinor = Math.sqrt(1 - (dEccentricity * dEccentricity)) * dMajor;

            rectOval.left = 0 - (float)(dMinor / 2);
            rectOval.right = (float)(dMinor / 2);

            rectOval.left += (float)iSunXOffset;
            rectOval.right += (float)iSunXOffset;
            rectOval.top += (float)iSunYOffset;
            rectOval.bottom += (float)iSunYOffset;

            // now correct for Longitude of Perihelion for this orbit's path

            canvas.save();
                canvas.rotate((float)dLongPeri, (float)iSunXOffset, (float)iSunYOffset);
                canvas.drawOval(rectOval, paint);
            canvas.restore();
            }

        int iBitmapHeight = bmp.getHeight();

        canvas.drawBitmap(bmp, fX - (iBitmapHeight / 2), fY - (iBitmapHeight / 2), null);

        // draw planet label

        myPaint.setColor(Color.WHITE);
        paint.setTextSize(30);
        canvas.drawText(sName, fX+20, fY-20, paint);
}

上述方法调用了另外两个方法，提供E（平均近点角）和r的值，r是行星所在向量的末端长度。 getE:

public double getE (double dTime, double dEccentricity)
    {
    // we are passed the degree count in degrees (duh) 
    // and the eccentricity value
    // the method returns E

    double dM1, dD, dE0, dE = 0; // return value E = the mean anomaly
    double dM; // local value of M in radians

    dM = Math.toRadians (dTime);

    int iSign = 1;

    if (dM > 0) iSign = 1; else iSign = -1;

    dM = Math.abs(dM) / (2 * Math.PI); // Meeus, p 206, line 110
    dM = (dM - (long)dM) * (2 * Math.PI) * iSign; // line 120
    if (dM < 0)
        dM = dM + (2 * Math.PI); // line 130
    iSign = 1;
    if (dM > Math.PI) iSign = -1; // line 150
    if (dM > Math.PI) dM = 2 * Math.PI - dM; // line 160

    dE0 = Math.PI / 2; // line 170
    dD = Math.PI / 4; // line 170

    for (int i = 0; i < 33; i++) // line 180 
        {
        dM1 = dE0 - dEccentricity * Math.sin(dE0); // line 190
        dE0 = dE0 + dD * Math.signum((float)(dM - dM1));
        dD = dD / 2; 
        }

    dE = dE0 * iSign;

    return dE;
    }

getRfromE:

public double getRfromE (double dE, double dEccentricty, double dScalar)
{
    return Math.min(getWidth(), getHeight()) / 2 * dScalar * (1 - (dEccentricty * Math.cos(dE)));
}

- Barry

很好的回答，如果这段代码是从书中复制/粘贴而来的，我们应该清楚地指出，以便给予那些赋予思想生命的作者应有的荣誉。 - Eric Leschinski

1

看看nMod，这是一个用C++编写并使用OpenGL的n体建模工具包。它有一个相当完整的太阳系模型，并且应该很容易修改。此外，他还有一个关于n体模拟的相当不错的维基百科。创建这个工具包的同一个人还在制作一个名为Moody的新程序，但似乎进展不如预期。

此外，如果您要进行超过几个对象的n体模拟，您真的应该看看快速多极子方法（也称为快速多极算法）。它可以将计算次数从O(N^2)降低到O(N)，以真正加速您的模拟。根据本文作者的说法，它也是20世纪最成功的十大算法之一之一。

- Justin Peel

1

看起来很难，需要对物理学有深厚的知识，但实际上它非常简单，你只需要知道两个公式和向量的基本理解：

质量为m1和m2，距离为d的行星1和行星2之间的引力（或重力）：Fg=G*m1*m2/d^2; Fg=m*a。G是一个常数，通过替换随机值找到它，使得加速度"a"不会太小也不会太大，大约为"0.01"或"0.1"。

如果你拥有作用于当前行星的总矢量力，你可以找到瞬时加速度a=(总力)/(当前行星的质量)。如果你有当前加速度、当前速度和当前位置，你可以找到新的速度和新的位置。

如果你想看到真实的效果，你可以使用以下超级简单的算法（伪代码）：

int n; // # of planets
Vector2D planetPosition[n]; 
Vector2D planetVelocity[n]; // initially set by (0, 0)
double planetMass[n];

while (true){
    for (int i = 0; i < n; i++){
        Vector2D totalForce = (0, 0); // acting on planet i
        for (int j = 0; j < n; j++){
            if (j == i)
                continue; // force between some planet and itself is 0
            Fg = G * planetMass[i] * planetMass[j] / distance(i, j) ^ 2;
        // Fg is a scalar value representing magnitude of force acting
        // between planet[i] and planet[j]
        // vectorFg is a vector form of force Fg
        // (planetPosition[j] - planetPosition[i]) is a vector value
        // (planetPosition[j]-planetPosition[i])/(planetPosition[j]-plantetPosition[i]).magnitude() is a
        // unit vector with direction from planet[i] to planet[j]
            vectorFg = Fg * (planetPosition[j] - planetPosition[i]) / 
                  (planetPosition[j] - planetPosition[i]).magnitude();
            totalForce += vectorFg;
        }
        Vector2D acceleration = totalForce / planetMass[i];
        planetVelocity[i] += acceleration;
    }

    // it is important to separate two for's, if you want to know why ask in the comments
    for (int i = 0; i < n; i++)
        planetPosition[i] += planetVelocity[i];

    sleep 17 ms;
    draw planets;
}

- rint

你能修复你的代码吗，特别是第vectorFg = magnitude is Fg, towards planet[j]行。我不确定它是否应该是一条注释，但因为它，我无法理解你的意思。 - kacpr

@kacpr 添加了全面的注释，希望能有所帮助！顺便说一下，这个算法没有考虑碰撞，如果你想要碰撞，需要添加一点代码。 - rint

你建议使用哪些行星质量和初始位置值作为代码的起点？我之前用过类似于你的代码，但当我输入一些以国际单位制给出的真实世界数据时，它就崩溃了。我在这里看到类似的情况发生，但我知道以千克为单位给出的太阳或地球的质量或以米为单位给出的距离是过度的。;) - kacpr

@kacpr 行星的距离和初始位置以像素为单位，大约相隔50到1000个像素，因此您可以在屏幕上看到它们。最好不要手动设置质量，可以根据半径在2D中使用行星的面积（或体积），例如：r = 50像素，质量= 2500 * pi。您唯一需要找出的是 G-万有引力常数，只需尝试将其设置为0.00001，然后为0.001，然后为10或1000，查看哪个值最合适。 - rint

0

如果你正在模拟物理现象，我强烈推荐使用Box2D。它是一个很棒的物理模拟器，能够大幅减少你在物理模拟中所需的样板代码。

- Joseph Salisbury

Box2D没有N体求解器或类似的东西。它假设一个恒定的重力势（至少截至2013年11月）。 - Domi

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Stefano Borini · Accepted Answer

这里包含了所有这里以及Jean Meeus所写的所有内容。

_{(来源: willbell.com)}