为纪念我的老师，新切尔卡斯克理工学院物理与数学学院第一任院长，亚历山大·尼古拉耶维奇·卡贝尔科夫（Alexander Nikolayevich Kabelkov）理论力学系主任

引言

八月，夏天快要结束了。人们疯狂地奔赴大海，是的，季节本身就不足为奇了。同时，在哈勃（Habr）上，伪科学开花并散发出强烈的色彩。如果我们谈论“建模...”这一主题，那么我们将把商业与快乐结合在一起-我们将继续承诺的循环，并用这种非常伪科学来与现代青年的好奇心作斗争。

但是这个问题并不是闲着吗？自从学年以来，我们就一直认为我们离外太空最近的卫星-月球绕地球运转了29.5天，特别是没有涉及到相关细节。实际上，我们的邻居是一种并且在某种程度上是独特的天文物体，它在地球上的移动并不像我的一些来自国外的同事所希望的那样简单。

因此，撇开争论，我们将在能力范围内从不同角度尝试考虑这一无条件的美丽，有趣和非常具有启发性的任务。

1.引力定律以及我们可以从中得出什么结论

引力定律是由艾萨克·牛顿爵士在17世纪下半叶发现的，它表明月球被吸引到地球（并且地球到月球！）的力沿着连接正在考虑的天体中心且大小相等的直线

， ，

$F_ {1,2} = G \，\ frac {m_1 \，m_2} {r_ {1,2} ^ 2}$

其中m ₁ ，m ₂分别是月球和地球的质量； G = 6.67e-11 m ³ /（kg * s ² ）-重力常数； r _1,2是月球中心与地球之间的距离。如果仅考虑这种力，那么，在解决了作为地球卫星的月亮运动的问题，并学会了在恒星背景下计算月亮在天空中的位置后，我们很快就会通过直接测量月球的赤道坐标而确信，在我们的音乐学院里，并非所有事物都如此光滑我想。而且这里的要点不是万有引力定律（在天体力学发展的早期，这种想法经常被表达），而是对月球运动对其他物体的不加激怒。哪一个我们看着天空，凝视着我们，紧紧地靠在鼻子下面的重1.99e30公斤等离子球-太阳下。月亮被太阳吸引了吗？就像绝对力相等

， ，

$F_ {1,3} = G \，\ frac {m_1 \，m_3} {r_ {1,3} ^ 2}$

m ₃是太阳的质量； r _1.3是从月亮到太阳的距离。将此功能与上一个功能进行比较。

， ， ， ， ， 左 （ 对 ）

$\ frac {F_ {1,3}} {F_ {1,2}} = \ frac {G \，\ frac {m_1 \，m_3} {r_ {1,3} ^ 2}} {G \，\ frac {m_1 \，m_2} {r_ {1,2} ^ 2}} = \ frac {m_3} {m_2} \，\左（\ frac {r_ {1,2}} {r_ {1,3}} \对）^ 2$

让我们采取这样的位置，在这种位置上，月亮对太阳的吸引力将最小：所有三个身体都在一条直线上，地球位于月亮和太阳之间。在这种情况下，我们的公式将采用以下形式：

， 左 （ ）

$\ frac {F_ {1,3}} {F_ {1,2}} = \ frac {m_3} {m_2} \，\左（\ frac {\ rho} {a + \ rho} \ right）^ 2$

在哪里

$\ rho = 3,844 \ cdot 10 ^ {8}$ ，m-地球到月球的平均距离；

$a = 1,496 \ cdot10 ^ {11}$ ，m-地球到太阳的平均距离。我们将实际参数代入该公式

$\ frac {F_ {1,3}} {F_ {1,2}} = \ frac {1.99 \ cdot 10 ^ {30}} {5.98 \ cdot10 ^ {24}} \，\ left（\ frac {3.844 \ cdot10 ^ {8}} {1.496 \ cdot10 ^ {11} + 3.844 \ cdot10 ^ {8}} \ right）^ 2 = $ 2.1$

这是号码！事实证明，月球被太阳吸引的力是其对地球吸引力的两倍以上。

这样的扰动不再被忽略，它将肯定影响月球的最终轨迹。让我们更进一步，考虑到地球轨道是半径为a的圆形轨道的假设，我们找到了地球周围各点的几何位置，其中任何物体对地球的吸引力等于其对太阳的吸引力。这将是一个半径为球形的球体

$R = \ frac {a \，\ sqrt {\ gamma}} {1-\ gamma}$

沿着连接地球和太阳的线向与太阳方向相反的一侧偏移一定距离

$l = R \，\ sqrt {\ gamma}$

在哪里

$\伽马= m_2 / m_3$ -地球质量与太阳质量之比。替换参数的数值，我们得出该区域的实际尺寸：R = 259300公里，l = 450公里。该球体称为地球相对于太阳的重力球体。

已知的月球轨道位于该区域之外。也就是说，在轨迹的任何一点上，月亮从太阳一侧比从地球一侧受到的吸引力要大得多。

2.卫星还是行星？引力范围

这些信息经常引起人们的争论，即月球不是地球的卫星，而是太阳系的独立行星，其轨道被附近地球的吸引力所扰动。

让我们使用拉普拉斯（P. Laplace）提出的标准，评估太阳对月球相对于地球的轨迹的扰动，以及地球对月球相对于太阳的轨迹的扰动。考虑三个物体：太阳（S），地球（E）和月球（M）。
我们假设地球相对于太阳的轨道和月亮相对于地球的轨道是圆形的。

考虑月球在地心惯性参考系中的运动。月球在日心参考系统中的绝对加速度取决于作用在月球上的重力，它等于：

$\ vec a_1 = \ vec a_1 ^ {（3）} + \ vec a_1 ^ {（2）} = \ frac {1} {m_1} \，\ vec F_ {1,3} + \ frac {1} { m_1} \，\ vec F_ {1,2}$

另一方面，根据科里奥利定理，月亮的绝对加速度

$\ vec a_1 = \ vec a_2 + \ vec a_ {1,2}$

在哪里

$\ vec a_2$ -便携式加速度等于地球相对于太阳的加速度；

$\ vec a_ {1,2}$ -月球相对于地球的加速度。这里不会有科里奥利加速度-我们选择的坐标系会逐渐移动。从这里我们得到月球相对于地球的加速度

$\ vec a_ {1,2} = \ frac {1} {m_1} \，\ vec F_ {1,3} + \ frac {1} {m_1} \，\ vec F_ {1,2}-\ vec a_2$

该加速度的一部分等于

$\ vec a_1 ^ {（（2）} = \ frac {1} {m_1} \，\ vec F_ {1,2}$ 由于月亮对地球的吸引力，并体现了其不受干扰的地心运动。其余的

$\ Delta \ vec a_ {1,3} = \ frac {1} {m_1} \，\ vec F_ {1,3}-\ vec a_2$

由太阳干扰引起的月球加速。

如果我们考虑日球在日心惯性参考系中的运动，那么一切都会简单得多，加速度

$\ vec a_1 ^ {（3）} = \ frac {1} {m_1} \，\ vec F_ {1,3}$ 表征了月亮不受干扰的日心运动和加速度

$\ Delta \ vec a_ {1,2} = \ frac {1} {m_1} \，\ vec F_ {1,2}$ -从地球一侧对该运动的干扰。

给定当前时代存在的地球和月球轨道的参数，不等式

$\ frac {| \ Delta \ vec a_ {1,3} |} {| \ vec a_ {1} ^ {（2）} |} <\ frac {| \ Delta \ vec a_ {1,2} |} {| \ vec a_ {1} ^ {（3）} |} \ Quad \ Quad（1）$

可以通过直接计算进行检查，但是我将参考资料来源，以免不必要地弄乱文章。

不平等（1）是什么意思？是的，相对而言，太阳对月球的扰动（非常明显）的影响要小于月球对地球的吸引力。相反，地球对月球的地心轨迹的愤慨对其运动的性质具有决定性的影响。在这种情况下，地球引力的影响更为显着，这意味着月亮“右”属于地球，是它的卫星。

另一件事很有趣-将不等式（1）转换为方程式，您可以找到月球（以及任何其他物体）的扰动对地球和太阳相同的点的几何位置。不幸的是，这并不像重力球那样简单。计算表明，该表面由疯狂阶方程描述，但接近于旋转椭圆体。我们要做的没有不必要麻烦的就是评估该表面相对于地球中心的整体尺寸。数值求解方程

$\ frac {| \ Delta \ vec a_ {1,3} |} {| \ vec a_ {1} ^ {（2）} |} = \ frac {| \ Delta \ vec a_ {1,2} |} {| \ vec a_ {1} ^ {（3）} |} \ Quad \ Quad（2）$

相对于在足够多的点上从地球中心到所需表面的距离，我们通过黄道面获得所需表面的横截面

为了清楚起见，此处显示了我们在上面相对于太阳发现的月球地心轨道和地球重力范围。从图中可以看出，影响范围或地球相对于太阳的重力作用范围是相对于X轴的旋转表面，沿连接地球和太阳的直线（沿日食轴）展平。月亮的轨道在这个假想的表面深处。

为了进行实际计算，可以方便地用一个球体近似该表面，该球体的中心位于地球中心，半径等于

$r =一个\，\左（\ frac {m} {M} \ right）^ {\ frac {2} {5}} \ quad \ quad（3）$

其中m是较小天体的质量； M是较大物体的质量，较小物体在其重力场中移动； a是身体中心之间的距离。就我们而言

$r =一个\，\左（\ frac {m_2} {m_3} \右）^ {\ frac {2} {5}} = 1.496 \ cdot10 ^ {11} \，\左（\ frac {5.98 \ cdot10 ^ {24}} {1.99 \ cdot10 ^ {30}} \ right）^ {\ frac {2} {5}} = 925000，\，km$

这个未完成的百万公里是理论极限，地球上的老太太的力量不会超出这个极限，她对天文物体轨迹的影响是如此之小，可以忽略不计。这意味着，在距离地球3840万公里的圆形轨道上发射月球（某些语言学家会这样做）会失败，从物理上讲是不可能的。

为了比较，该球体在图中用蓝色虚线表示。在评估计算中，假设该球体内部的物体将完全从地球一侧经历引力。如果身体在此球体之外，则我们认为该身体在太阳的引力场中移动。在实际的航空航天中，共轭圆锥形截面的方法是已知的，它允许使用二体问题的解近似地计算航天器的轨迹。此外，设备克服的所有空间都分为相似的影响范围。

例如，现在很明显，为了具有理论上的能力使操纵进入近月轨道，航天器必须落在月球相对于地球的作用范围内。它的半径很容易通过公式（3）计算出来，为6.6万公里。

因此，月球可以正确地视为地球的卫星。但是，由于太阳引力场的重大影响，它不会在中心引力场中移动，这意味着它的轨迹不是圆锥形截面。

3.经典公式中的三体问题

因此，我们将在一般情况下考虑模型问题，在天体力学中称为三体问题。让我们考虑任意质量的三个物体，它们任意地位于空间中，并且仅在相互引力的作用下移动

我们认为身体是物质点。物体的位置将在任意基础上计算，与惯性参考系统Oxyz相关联。每个物体的位置分别由半径矢量设置

$\ vec r_1$ ，

$\ vec r_2$ 和

$\ vec r_3$ 。此外，根据点动力学的第三公理（牛顿第三定律），来自其他两个物体侧面的引力也作用在每个物体上。

$\ vec F_ {i，j} =-\ vec F_ {j，i} \ quad \ quad（4）$

我们以矢量形式编写每个点的运动微分方程

$\开始{align}＆m_1 \，\ frac {d ^ 2 \ vec r_1} {dt ^ 2} = \ vec F_ {1,2} + \ vec F_ {1,3} \\＆m_2 \，\ frac {d ^ 2 \ vec r_2} {dt ^ 2} = \ vec F_ {2,1} + \ vec F_ {2,3} \\＆m_3 \，\ frac {d ^ 2 \ vec r_3} {dt ^ 2} = \ vec F_ {3,1} + \ vec F_ {3,2} \ end {align}$

或者，考虑到（4）

$\开始{align}＆m_1 \，\ frac {d ^ 2 \ vec r_1} {dt ^ 2} = \ vec F_ {1,2} + \ vec F_ {1,3} \\＆m_2 \，\ frac {d ^ 2 \ vec r_2} {dt ^ 2} =-\ vec F_ {1,2} + \ vec F_ {2,3} \\＆m_3 \，\ frac {d ^ 2 \ vec r_3} { dt ^ 2} =-\ vec F_ {1,3}-\ vec F_ {2,3} \ end {align}$

根据万有引力定律，相互作用力沿矢量定向

$\开始{align}＆\ vec r_ {1,2} = \ vec r_2-\ vec r_1 \\＆\ vec r_ {1,3} = \ vec r_3-\ vec r_1 \\＆\ vec r_ {2 ，3} = \ vec r_3-\ vec r_2 \\ \结束{align}$

沿着每个向量，我们发布一个对应的单位向量

$\ vec e_ {i，j} = \ frac {1} {r_ {i，j}} \，\ vec r_ {i，j}$

然后每个重力由下式计算

$\ vec F_ {i，j} = G \，\ frac {m_i \，m_j} {r_ {i，j} ^ 2} \，\ vec e_ {i，j}$

鉴于所有这些，运动方程组采用以下形式：

$\开始{align}＆\ frac {d ^ 2 \ vec r_1} {dt ^ 2} = \ frac {G \，m_2} {r_ {1,2} ^ 3} \，\ vec r_ {1,2 } + \ frac {G \，m_3} {r_ {1,3} ^ 3} \，\ vec r_ {1,3} \\＆\ frac {d ^ 2 \ vec r_2} {dt ^ 2} =- \ frac {G \，m_1} {r_ {1,2} ^ 3} \，\ vec r_ {1,2} + \ frac {G \，m_3} {r_ {2,3} ^ 3} \，\ vec r_ {2,3} \\＆\ frac {d ^ 2 \ vec r_3} {dt ^ 2} =-\ frac {G \，m_1} {r_ {1,3} ^ 3} \，\ vec r_ {1,3}-\ frac {G \，m_2} {r_ {2,3} ^ 3} \，\ vec r_ {2,3} \ end {align}$

介绍天体力学接受的符号

$\ mu_i = G \，m_i$

是引力中心的引力参数。然后，运动方程将采用最终的矢量形式

$\开始{align}＆\ frac {d ^ 2 \ vec r_1} {dt ^ 2} = \ frac {\ mu_2} {r_ {1,2} ^ 3} \，\ vec r_ {1,2} + \ frac {\ mu_3} {r_ {1,3} ^ 3} \，\ vec r_ {1,3} \\＆\ frac {d ^ 2 \ vec r_2} {dt ^ 2} =-\ frac {\ mu_1} {r_ {1,2} ^ 3} \，\ vec r_ {1,2} + \ frac {\ mu_3} {r_ {2,3} ^ 3} \，\ vec r_ {2,3} \ \＆\ frac {d ^ 2 \ vec r_3} {dt ^ 2} =-\ frac {\ mu_1} {r_ {1,3} ^ 3} \，\ vec r_ {1,3}-\ frac {\ mu_2} {r_ {2,3} ^ 3} \，\ vec r_ {2,3} \ end {align}$

4.方程与无量纲变量的比例分配

数学建模中一种相当流行的技术是将微分方程和其他关系描述为无因次相坐标和无因次时间的关系。其他参数也被标准化。尽管使用了数值模型，这使我们可以考虑，但以相当普遍的形式考虑了一整套典型问题。我保留了在每个要解决的问题中这样做有多合理的问题，但是我同意在这种情况下这种方法是相当公平的。

因此，我们介绍了一些具有引力参数的抽象天体

$\亩$ ，以使卫星在具有长半轴的椭圆形轨道上的自转周期

$一$ 周围等于

$T$ 。所有这些量，根据力学定律，都通过以下关系关联

$T = 2 \，\ pi \，\左（\ frac {a ^ 3} {\ mu} \ right）^ {\ frac {1} {2}}$

我们介绍参数的替换。对于我们系统中各点的位置

$\ vec r_i =一个\，\ vec \ xi_i$

在哪里

$\ vec \ xi_i$ -第i个点的无量纲半径矢量；
物体的重力参数

$\ mu_i = \ varkappa_i \，\ mu$

在哪里

$\ varkappa_i$ -第i个点的无量纲引力参数；
时间

$t = T \，\ tau$

在哪里

$\ tau$ -无量纲的时间。

现在，我们通过这些无量纲参数重新计算系统的加速点。我们应用直接两倍时间微分。为了速度

$\ vec v_i = \ frac {d \ vec r_i} {dt} = a \，\ frac {d \ vec \ xi_i} {dt} = \ frac {a} {T} \，\ frac {d \ vec \ xi_i} {d \ tau} = \ frac {1} {2 \，\ pi} \，\ sqrt {\ frac {\ mu} {a}} \，\ frac {d \ vec \ xi_i} {d \ tau }。$

为了加速

$\ vec a_i = \ frac {d \ vec v_i} {dt} = \ frac {1} {2 \，\ pi} \，\ sqrt {\ frac {\ mu} {a}} \，\ frac {1 } {dt} \左（\ frac {d \ vec \ xi_i} {d \ tau} \右）= \ frac {1} {4 \，\ pi ^ 2} \，\ frac {\ mu} {a ^ 2} \，\ frac {d ^ 2 \ vec \ xi_i} {d \ tau ^ 2}$

当将获得的关系代入运动方程式时，一切都会优雅地折叠成漂亮的方程式：

$\ begin {align}＆\ frac {d ^ 2 \ vec \ xi_1} {d \ tau ^ 2} = 4 \，\ pi ^ 2 \，\ varkappa_2 \，\ frac {\ vec \ xi_2-\ vec \ xi_1} {| \ vec \ xi_2-\ vec \ xi_1 | ^ 3} + 4 \，\ pi ^ 2 \，\ varkappa_3 \，\ frac {\ vec \ xi_3-\ vec \ xi_1} {| \ vec \ xi_3 -\ vec \ xi_1 | ^ 3} \\＆\ frac {d ^ 2 \ vec \ xi_2} {d \ tau ^ 2} = -4 \，\ pi ^ 2 \，\ varkappa_1 \，\ frac {\ vec \ xi_2-\ vec \ xi_1} {| \ vec \ xi_2-\ vec \ xi_1 | ^ 3} + 4 \，\ pi ^ 2 \，\ varkappa_3 \，\ frac {\ vec \ xi_3-\ vec \ xi_2} {| \ vec \ xi_3-\ vec \ xi_2 | ^ 3} \ quad \ quad（5）\\＆\ frac {d ^ 2 \ vec \ xi_3} {d \ tau ^ 2} = -4 \，\ pi ^ 2 \，\ varkappa_1 \，\ frac {\ vec \ xi_3-\ vec \ xi_1} {| \ vec \ xi_3-\ vec \ xi_1 | ^ 3}-4 \，\ pi ^ 2 \，\ varkappa_2 \， \ frac {\ vec \ xi_3-\ vec \ xi_2} {| \ vec \ xi_3-\ vec \ xi_2 | ^ 3} \ end {align}$

该方程组仍被认为在解析函数中不可积分。为什么要考虑而不是呢？由于复变函数理论的成功导致了三体问题的通用解在1912年出现的事实-卡尔·宗德曼（Karl Zundman）找到了一种算法，用于针对复杂参数找到无限级数的系数，从理论上讲，这是三体问题的通用解。但是...为了使Sundman系列以其所需的精度在实际计算中的应用，它需要获得这些系列的成员数量，以至于即使在今天，这项任务也远远超出了计算机的功能。

因此，数值积分是分析方程式（5）解的唯一方法

5.初始条件的计算：我们获得了初始数据

正如我之前所写，在开始数值积分之前，您应该注意计算要解决的问题的初始条件。在所考虑的问题中，搜索初始条件变成一个独立的子问题，因为系统（5）给出了九个二阶标量方程，当传递给标准柯西形式时，系统的阶数增加了2倍。也就是说，我们需要计算多达18个参数-系统中所有点的初始位置和初始速度的分量。我们从哪里获得有关我们感兴趣的天体位置的数据？我们生活在一个人在月球上行走的世界-自然，人类必须掌握有关月球如何移动以及它位于何处的信息。

也就是说，您说，伙计，您是在为我们提供从架子上取下厚厚的天文书籍，将它们吹掉的灰尘……别猜！我建议去那些实际上在月球上行走的人，即加利福尼亚州帕萨迪纳的喷气推进实验室去NASA。此处-JPL Horizonts Web界面。

在这里，花了一些时间研究接口之后，我们将获得所需的所有数据。选择一个日期，例如，是的，我们不在乎，但是将其设为2018年7月27日UT 20:21。就在这时，观测到月全食。该计划将为我们提供巨大的基础

2018年7月27日20:21月球星历的完整结论（起源于地球中心）

******************************************************************************* Revised: Jul 31, 2013 Moon / (Earth) 301 GEOPHYSICAL DATA (updated 2018-Aug-13): Vol. Mean Radius, km = 1737.53+-0.03 Mass, x10^22 kg = 7.349 Radius (gravity), km = 1738.0 Surface emissivity = 0.92 Radius (IAU), km = 1737.4 GM, km^3/s^2 = 4902.800066 Density, g/cm^3 = 3.3437 GM 1-sigma, km^3/s^2 = +-0.0001 V(1,0) = +0.21 Surface accel., m/s^2 = 1.62 Earth/Moon mass ratio = 81.3005690769 Farside crust. thick. = ~80 - 90 km Mean crustal density = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km Heat flow, Apollo 15 = 3.1+-.6 mW/m^2 k2 = 0.024059 Heat flow, Apollo 17 = 2.2+-.5 mW/m^2 Rot. Rate, rad/s = 0.0000026617 Geometric Albedo = 0.12 Mean angular diameter = 31'05.2" Orbit period = 27.321582 d Obliquity to orbit = 6.67 deg Eccentricity = 0.05490 Semi-major axis, a = 384400 km Inclination = 5.145 deg Mean motion, rad/s = 2.6616995x10^-6 Nodal period = 6798.38 d Apsidal period = 3231.50 d Mom. of inertia C/MR^2= 0.393142 beta (CA/B), x10^-4 = 6.310213 gamma (BA/C), x10^-4 = 2.277317 Perihelion Aphelion Mean Solar Constant (W/m^2) 1414+-7 1323+-7 1368+-7 Maximum Planetary IR (W/m^2) 1314 1226 1268 Minimum Planetary IR (W/m^2) 5.2 5.2 5.2 ******************************************************************************* ******************************************************************************* Ephemeris / WWW_USER Wed Aug 15 20:45:05 2018 Pasadena, USA / Horizons ******************************************************************************* Target body name: Moon (301) {source: DE431mx} Center body name: Earth (399) {source: DE431mx} Center-site name: BODY CENTER ******************************************************************************* Start time : AD 2018-Jul-27 20:21:00.0003 TDB Stop time : AD 2018-Jul-28 20:21:00.0003 TDB Step-size : 0 steps ******************************************************************************* Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)} Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)} Center radii : 6378.1 x 6378.1 x 6356.8 km {Equator, meridian, pole} Output units : AU-D Output type : GEOMETRIC cartesian states Output format : 3 (position, velocity, LT, range, range-rate) Reference frame : ICRF/J2000.0 Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ******************************************************************************* JDTDB XYZ VX VY VZ LT RG RR ******************************************************************************* $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 1.537109094089627E-03 Y =-2.237488447258137E-03 Z = 5.112037386426180E-06 VX= 4.593816208618667E-04 VY= 3.187527302531735E-04 VZ=-5.183707711777675E-05 LT= 1.567825598846416E-05 RG= 2.714605874095336E-03 RR=-2.707898607099066E-06 $$EOE ******************************************************************************* Coordinate system description: Ecliptic and Mean Equinox of Reference Epoch Reference epoch: J2000.0 XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76) X-axis : out along ascending node of instantaneous plane of the Earth's orbit and the Earth's mean equator at the reference epoch Z-axis : perpendicular to the xy-plane in the directional (+ or -) sense of Earth's north pole at the reference epoch. Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]: JDTDB Julian Day Number, Barycentric Dynamical Time X X-component of position vector (au) Y Y-component of position vector (au) Z Z-component of position vector (au) VX X-component of velocity vector (au/day) VY Y-component of velocity vector (au/day) VZ Z-component of velocity vector (au/day) LT One-way down-leg Newtonian light-time (day) RG Range; distance from coordinate center (au) RR Range-rate; radial velocity wrt coord. center (au/day) Geometric states/elements have no aberrations applied. Computations by ... Solar System Dynamics Group, Horizons On-Line Ephemeris System 4800 Oak Grove Drive, Jet Propulsion Laboratory Pasadena, CA 91109 USA Information: http://ssd.jpl.nasa.gov/ Connect : telnet://ssd.jpl.nasa.gov:6775 (via browser) http://ssd.jpl.nasa.gov/?horizons telnet ssd.jpl.nasa.gov 6775 (via command-line) Author : Jon.D.Giorgini@jpl.nasa.gov *******************************************************************************

Brrr，这是什么？如果没有恐慌，对于在学校里精通天文学，力学和数学的人来说，没有什么可担心的。因此，最重要的是月球最终寻找的坐标和速度分量。

 $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 1.537109094089627E-03 Y =-2.237488447258137E-03 Z = 5.112037386426180E-06 VX= 4.593816208618667E-04 VY= 3.187527302531735E-04 VZ=-5.183707711777675E-05 LT= 1.567825598846416E-05 RG= 2.714605874095336E-03 RR=-2.707898607099066E-06 $$EOE

是的，是的，他们是笛卡尔！如果仔细阅读整个鞋底布，我们会发现该坐标系的原点与地球中心重合。 XY平面位于J2000时代的地球轨道平面（黄道平面）中。 X轴沿着地球的赤道平面与黄道点在黄道的相交线指向。 Z轴的方向是垂直于黄道平面的地球北极。好吧，Y轴将所有这些快乐补充到正确的三个向量上。默认情况下，坐标单位：天文单位（来自NASA的数字也给出了以千米为单位的自治单位的大小）。速度单位：每天的天文单位，假定天为86400秒。全馅！

我们可以获得有关地球的类似信息。

地球星历的完整结论于07/27/2018 20:21（原点位于太阳系质心）

 ******************************************************************************* Revised: Jul 31, 2013 Earth 399 GEOPHYSICAL PROPERTIES (revised Aug 13, 2018): Vol. Mean Radius (km) = 6371.01+-0.02 Mass x10^24 (kg)= 5.97219+-0.0006 Equ. radius, km = 6378.137 Mass layers: Polar axis, km = 6356.752 Atmos = 5.1 x 10^18 kg Flattening = 1/298.257223563 oceans = 1.4 x 10^21 kg Density, g/cm^3 = 5.51 crust = 2.6 x 10^22 kg J2 (IERS 2010) = 0.00108262545 mantle = 4.043 x 10^24 kg g_p, m/s^2 (polar) = 9.8321863685 outer core = 1.835 x 10^24 kg g_e, m/s^2 (equatorial) = 9.7803267715 inner core = 9.675 x 10^22 kg g_o, m/s^2 = 9.82022 Fluid core rad = 3480 km GM, km^3/s^2 = 398600.435436 Inner core rad = 1215 km GM 1-sigma, km^3/s^2 = 0.0014 Escape velocity = 11.186 km/s Rot. Rate (rad/s) = 0.00007292115 Surface Area: Mean sidereal day, hr = 23.9344695944 land = 1.48 x 10^8 km Mean solar day 2000.0, s = 86400.002 sea = 3.62 x 10^8 km Mean solar day 1820.0, s = 86400.0 Moment of inertia = 0.3308 Love no., k2 = 0.299 Mean Temperature, K = 270 Atm. pressure = 1.0 bar Vis. mag. V(1,0) = -3.86 Volume, km^3 = 1.08321 x 10^12 Geometric Albedo = 0.367 Magnetic moment = 0.61 gauss Rp^3 Solar Constant (W/m^2) = 1367.6 (mean), 1414 (perihelion), 1322 (aphelion) ORBIT CHARACTERISTICS: Obliquity to orbit, deg = 23.4392911 Sidereal orb period = 1.0000174 y Orbital speed, km/s = 29.79 Sidereal orb period = 365.25636 d Mean daily motion, deg/d = 0.9856474 Hill's sphere radius = 234.9 ******************************************************************************* ******************************************************************************* Ephemeris / WWW_USER Wed Aug 15 21:16:21 2018 Pasadena, USA / Horizons ******************************************************************************* Target body name: Earth (399) {source: DE431mx} Center body name: Solar System Barycenter (0) {source: DE431mx} Center-site name: BODY CENTER ******************************************************************************* Start time : AD 2018-Jul-27 20:21:00.0003 TDB Stop time : AD 2018-Jul-28 20:21:00.0003 TDB Step-size : 0 steps ******************************************************************************* Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)} Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)} Center radii : (undefined) Output units : AU-D Output type : GEOMETRIC cartesian states Output format : 3 (position, velocity, LT, range, range-rate) Reference frame : ICRF/J2000.0 Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ******************************************************************************* JDTDB XYZ VX VY VZ LT RG RR ******************************************************************************* $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 5.755663665315949E-01 Y =-8.298818915224488E-01 Z =-5.366994499016168E-05 VX= 1.388633512282171E-02 VY= 9.678934168415631E-03 VZ= 3.429889230737491E-07 LT= 5.832932117417083E-03 RG= 1.009940888883960E+00 RR=-3.947237246302148E-05 $$EOE ******************************************************************************* Coordinate system description: Ecliptic and Mean Equinox of Reference Epoch Reference epoch: J2000.0 XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76) X-axis : out along ascending node of instantaneous plane of the Earth's orbit and the Earth's mean equator at the reference epoch Z-axis : perpendicular to the xy-plane in the directional (+ or -) sense of Earth's north pole at the reference epoch. Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]: JDTDB Julian Day Number, Barycentric Dynamical Time X X-component of position vector (au) Y Y-component of position vector (au) Z Z-component of position vector (au) VX X-component of velocity vector (au/day) VY Y-component of velocity vector (au/day) VZ Z-component of velocity vector (au/day) LT One-way down-leg Newtonian light-time (day) RG Range; distance from coordinate center (au) RR Range-rate; radial velocity wrt coord. center (au/day) Geometric states/elements have no aberrations applied. Computations by ... Solar System Dynamics Group, Horizons On-Line Ephemeris System 4800 Oak Grove Drive, Jet Propulsion Laboratory Pasadena, CA 91109 USA Information: http://ssd.jpl.nasa.gov/ Connect : telnet://ssd.jpl.nasa.gov:6775 (via browser) http://ssd.jpl.nasa.gov/?horizons telnet ssd.jpl.nasa.gov 6775 (via command-line) Author : Jon.D.Giorgini@jpl.nasa.gov *******************************************************************************

在此，将太阳系的重心（质心）选择为原点。我们感兴趣的数据

 $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 5.755663665315949E-01 Y =-8.298818915224488E-01 Z =-5.366994499016168E-05 VX= 1.388633512282171E-02 VY= 9.678934168415631E-03 VZ= 3.429889230737491E-07 LT= 5.832932117417083E-03 RG= 1.009940888883960E+00 RR=-3.947237246302148E-05 $$EOE

对于月球，我们需要相对于太阳系重心的坐标和速度，我们可以对其进行计算，也可以要求NASA提供此类数据

月球星历的完整结论于07/27/2018 20:21（太阳系质心的起源）

 ******************************************************************************* Revised: Jul 31, 2013 Moon / (Earth) 301 GEOPHYSICAL DATA (updated 2018-Aug-13): Vol. Mean Radius, km = 1737.53+-0.03 Mass, x10^22 kg = 7.349 Radius (gravity), km = 1738.0 Surface emissivity = 0.92 Radius (IAU), km = 1737.4 GM, km^3/s^2 = 4902.800066 Density, g/cm^3 = 3.3437 GM 1-sigma, km^3/s^2 = +-0.0001 V(1,0) = +0.21 Surface accel., m/s^2 = 1.62 Earth/Moon mass ratio = 81.3005690769 Farside crust. thick. = ~80 - 90 km Mean crustal density = 2.97+-.07 g/cm^3 Nearside crust. thick.= 58+-8 km Heat flow, Apollo 15 = 3.1+-.6 mW/m^2 k2 = 0.024059 Heat flow, Apollo 17 = 2.2+-.5 mW/m^2 Rot. Rate, rad/s = 0.0000026617 Geometric Albedo = 0.12 Mean angular diameter = 31'05.2" Orbit period = 27.321582 d Obliquity to orbit = 6.67 deg Eccentricity = 0.05490 Semi-major axis, a = 384400 km Inclination = 5.145 deg Mean motion, rad/s = 2.6616995x10^-6 Nodal period = 6798.38 d Apsidal period = 3231.50 d Mom. of inertia C/MR^2= 0.393142 beta (CA/B), x10^-4 = 6.310213 gamma (BA/C), x10^-4 = 2.277317 Perihelion Aphelion Mean Solar Constant (W/m^2) 1414+-7 1323+-7 1368+-7 Maximum Planetary IR (W/m^2) 1314 1226 1268 Minimum Planetary IR (W/m^2) 5.2 5.2 5.2 ******************************************************************************* ******************************************************************************* Ephemeris / WWW_USER Wed Aug 15 21:19:24 2018 Pasadena, USA / Horizons ******************************************************************************* Target body name: Moon (301) {source: DE431mx} Center body name: Solar System Barycenter (0) {source: DE431mx} Center-site name: BODY CENTER ******************************************************************************* Start time : AD 2018-Jul-27 20:21:00.0003 TDB Stop time : AD 2018-Jul-28 20:21:00.0003 TDB Step-size : 0 steps ******************************************************************************* Center geodetic : 0.00000000,0.00000000,0.0000000 {E-lon(deg),Lat(deg),Alt(km)} Center cylindric: 0.00000000,0.00000000,0.0000000 {E-lon(deg),Dxy(km),Dz(km)} Center radii : (undefined) Output units : AU-D Output type : GEOMETRIC cartesian states Output format : 3 (position, velocity, LT, range, range-rate) Reference frame : ICRF/J2000.0 Coordinate systm: Ecliptic and Mean Equinox of Reference Epoch ******************************************************************************* JDTDB XYZ VX VY VZ LT RG RR ******************************************************************************* $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 5.771034756256845E-01 Y =-8.321193799697072E-01 Z =-4.855790760378579E-05 VX= 1.434571674368357E-02 VY= 9.997686898668805E-03 VZ=-5.149408819470315E-05 LT= 5.848610189172283E-03 RG= 1.012655462859054E+00 RR=-3.979984423450087E-05 $$EOE ******************************************************************************* Coordinate system description: Ecliptic and Mean Equinox of Reference Epoch Reference epoch: J2000.0 XY-plane: plane of the Earth's orbit at the reference epoch Note: obliquity of 84381.448 arcseconds wrt ICRF equator (IAU76) X-axis : out along ascending node of instantaneous plane of the Earth's orbit and the Earth's mean equator at the reference epoch Z-axis : perpendicular to the xy-plane in the directional (+ or -) sense of Earth's north pole at the reference epoch. Symbol meaning [1 au= 149597870.700 km, 1 day= 86400.0 s]: JDTDB Julian Day Number, Barycentric Dynamical Time X X-component of position vector (au) Y Y-component of position vector (au) Z Z-component of position vector (au) VX X-component of velocity vector (au/day) VY Y-component of velocity vector (au/day) VZ Z-component of velocity vector (au/day) LT One-way down-leg Newtonian light-time (day) RG Range; distance from coordinate center (au) RR Range-rate; radial velocity wrt coord. center (au/day) Geometric states/elements have no aberrations applied. Computations by ... Solar System Dynamics Group, Horizons On-Line Ephemeris System 4800 Oak Grove Drive, Jet Propulsion Laboratory Pasadena, CA 91109 USA Information: http://ssd.jpl.nasa.gov/ Connect : telnet://ssd.jpl.nasa.gov:6775 (via browser) http://ssd.jpl.nasa.gov/?horizons telnet ssd.jpl.nasa.gov 6775 (via command-line) Author : Jon.D.Giorgini@jpl.nasa.gov *******************************************************************************

 $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 5.771034756256845E-01 Y =-8.321193799697072E-01 Z =-4.855790760378579E-05 VX= 1.434571674368357E-02 VY= 9.997686898668805E-03 VZ=-5.149408819470315E-05 LT= 5.848610189172283E-03 RG= 1.012655462859054E+00 RR=-3.979984423450087E-05 $$EOE

太好了！现在，您需要使用文件稍微处理接收到的数据。

6.38只鹦鹉和一只鹦鹉翼

首先，我们将确定比例，因为我们的运动方程（5）以无量纲形式书写。 NASA本身提供的数据告诉我们，值得以一个天文单位作为坐标比例尺。因此，我们将以太阳为标准物体，将其标准化为其他物体的质量，以地球绕太阳公转的周期为时间尺度。

当然，所有这些都很好，但是我们没有为太阳设定初始条件。 “为什么？” 语言学家会问我。我会回答，太阳绝不会静止，而是围绕太阳系质心在其轨道中旋转。通过查看太阳的NASA数据可以看出这一点。

 $$SOE 2458327.347916670 = AD 2018-Jul-27 20:21:00.0003 TDB X = 6.520050993518213E+04 Y = 1.049687363172734E+06 Z =-1.304404963058507E+04 VX=-1.265326939350981E-02 VY= 5.853475278436883E-03 VZ= 3.136673455633667E-04 LT= 3.508397935601254E+00 RG= 1.051791240756026E+06 RR= 5.053500842402456E-03 $$EOE

查看RG参数，我们看到太阳围绕太阳系的重心旋转，并且在2018年7月27日，恒星的中心距离它有一百万公里。太阳半径，仅供参考-69.6万公里。也就是说，太阳系的重心距离恒星表面五百万公里。怎么了是的，因为与太阳相互作用的所有其他物体也会使它加速，当然主要是沉重的木星。因此，太阳也有自己的轨道。

当然，我们可以选择此数据作为初始条件，但不能-我们正在解决模型三体问题，并且不包括木星和其他字符。因此，在知道了地球和月球的位置和速度的情况下，有损于现实主义，我们重新计算了太阳的初始条件，从而使太阳-地球-月亮系统的质心在原点。对于我们机械系统的质心，方程

$（m_1 + m_2 + m_3）\，\ vec r_C = m_1 \，\ vec r_1 + m_2 \，\ vec r_2 + m_3 \，\ vec r_3$

我们将质心放在原点，也就是说，我们问

$\ vec r_C = 0$ 然后

$m_1 \，\ vec r_1 + m_2 \，\ vec r_2 + m_3 \，\ vec r_3 = 0$

从哪里来

$\开始{align}＆m_3 \，\ vec r_3 = -m_1 \，\ vec r_1-m_2 \，\ vec r_2 \\＆\ vec r_3 =-\ frac {m_1} {m_3} \ vec r_1-\ frac {m_2} {m_3} \，\ vec r_2 \结束{align}$

让我们通过选择继续到无量纲坐标和参数

$\ mu = \ mu_3$

$\ vec \ xi_3 =-\ varkappa_1 \ vec \ xi_1-\ varkappa_2 \ vec \ xi_2 \ quad \ quad（6）$

关于时间微分（6）并传递到无量纲时间，我们还获得了速度的关系

$\ vec u_3 =-\ varkappa_1 \，\ vec u_1-\ varkappa_2 \，\ vec u_2$

在哪里

$\ vec u_i = \ cfrac {d \ vec \ xi_i} {d \ tau}，\ forall i = \ overline {1,3}$

现在，我们将编写一个程序，该程序将在我们选择的“鹦鹉”中形成初始条件。我们会写些什么？当然在Python中！毕竟，如您所知，这是数学建模的最佳语言。

但是，如果我们摆脱讽刺，那么我们真的会为此尝试python，为什么不呢？我一定会在我的Github个人资料中提供指向所有代码的链接。

月球-地球-太阳系初始条件的计算。

 # #    # #   G = 6.67e-11 #   (, , ) m = [7.349e22, 5.792e24, 1.989e30] #     mu = [] print("  ") for i, mass in enumerate(m): mu.append(G * mass) print("mu[" + str(i) + "] = " + str(mu[i])) #      kappa = [] print("  ") for i, gp in enumerate(mu): kappa.append(gp / mu[2]) print("xi[" + str(i) + "] = " + str(kappa[i])) print("\n") #   a = 1.495978707e11 import math #   , c T = 2 * math.pi * a * math.sqrt(a / mu[2]) print("  T = " + str(T) + "\n") #  NASA   xL = 5.771034756256845E-01 yL = -8.321193799697072E-01 zL = -4.855790760378579E-05 import numpy as np xi_10 = np.array([xL, yL, zL]) print("  , ..: " + str(xi_10)) #  NASA   xE = 5.755663665315949E-01 yE = -8.298818915224488E-01 zE = -5.366994499016168E-05 xi_20 = np.array([xE, yE, zE]) print("  , ..: " + str(xi_20)) #    ,     -      xi_30 = - kappa[0] * xi_10 - kappa[1] * xi_20 print("  , ..: " + str(xi_30)) #       Td = 86400.0 u = math.sqrt(mu[2] / a) / 2 / math.pi print("\n") #    vxL = 1.434571674368357E-02 vyL = 9.997686898668805E-03 vzL = -5.149408819470315E-05 vL0 = np.array([vxL, vyL, vzL]) uL0 = np.array([0.0, 0.0, 0.0]) for i, v in enumerate(vL0): vL0[i] = v * a / Td uL0[i] = vL0[i] / u print("  , /: " + str(vL0)) print(" -//- : " + str(uL0)) #    vxE = 1.388633512282171E-02 vyE = 9.678934168415631E-03 vzE = 3.429889230737491E-07 vE0 = np.array([vxE, vyE, vzE]) uE0 = np.array([0.0, 0.0, 0.0]) for i, v in enumerate(vE0): vE0[i] = v * a / Td uE0[i] = vE0[i] / u print("  , /: " + str(vE0)) print(" -//- : " + str(uE0)) #    vS0 = - kappa[0] * vL0 - kappa[1] * vE0 uS0 = - kappa[0] * uL0 - kappa[1] * uE0 print("  , /: " + str(vS0)) print(" -//- : " + str(uS0))

程序排气

    mu[0] = 4901783000000.0 mu[1] = 386326400000000.0 mu[2] = 1.326663e+20    xi[0] = 3.6948215183509304e-08 xi[1] = 2.912016088486677e-06 xi[2] = 1.0   T = 31563683.35432583   , ..: [ 5.77103476e-01 -8.32119380e-01 -4.85579076e-05]   , ..: [ 5.75566367e-01 -8.29881892e-01 -5.36699450e-05]   , ..: [-1.69738146e-06 2.44737475e-06 1.58081871e-10]   , /: [24838.98933473 17310.56333294 -89.15979106] -//- : [ 5.24078311 3.65235907 -0.01881184]   , /: [2.40435899e+04 1.67586567e+04 5.93870516e-01] -//- : [5.07296163e+00 3.53591219e+00 1.25300854e-04]   , /: [-7.09330769e-02 -4.94410725e-02 1.56493465e-06] -//- : [-1.49661835e-05 -1.04315813e-05 3.30185861e-10]

7.运动方程式的积分和结果分析

实际上，积分本身或多或少地简化为准备SciPy方程组的标准过程：将ODE系统转换为Cauchy形式并调用相应的求解器函数。要将系统转换为柯西格式，我们记得

$\ vec u_i = \ frac {d \ vec \ xi_i} {d \ tau}，\ forall i = \ overline {1,3} \ quad \ quad（7）$

然后介绍系统的状态向量

$\ vec y = \左[\ vec \ xi_1，\ vec \ xi_2，\ vec \ xi_1，\ vec u_1，\ vec u_2，\ vec u_3 \右侧] ^ T$

我们将（7）和（5）简化为一个向量方程

$\ frac {d \ vec y} {d \ tau} = \ vec f（\ tau，\ vec y）\ quad \ quad（8）$

为了将（8）与现有的初始条件集成在一起，我们编写了一些非常少的代码

三体问题中运动方程的积分

 # #     # def calcAccels(xi): k = 4 * math.pi ** 2 xi12 = xi[1] - xi[0] xi13 = xi[2] - xi[0] xi23 = xi[2] - xi[1] s12 = math.sqrt(np.dot(xi12, xi12)) s13 = math.sqrt(np.dot(xi13, xi13)) s23 = math.sqrt(np.dot(xi23, xi23)) a1 = (k * kappa[1] / s12 ** 3) * xi12 + (k * kappa[2] / s13 ** 3) * xi13 a2 = -(k * kappa[0] / s12 ** 3) * xi12 + (k * kappa[2] / s23 ** 3) * xi23 a3 = -(k * kappa[0] / s13 ** 3) * xi13 - (k * kappa[1] / s23 ** 3) * xi23 return [a1, a2, a3] # #       # def f(t, y): n = 9 dydt = np.zeros((2 * n)) for i in range(0, n): dydt[i] = y[i + n] xi1 = np.array(y[0:3]) xi2 = np.array(y[3:6]) xi3 = np.array(y[6:9]) accels = calcAccels([xi1, xi2, xi3]) i = n for accel in accels: for a in accel: dydt[i] = a i = i + 1 return dydt #     y0 = [xi_10[0], xi_10[1], xi_10[2], xi_20[0], xi_20[1], xi_20[2], xi_30[0], xi_30[1], xi_30[2], uL0[0], uL0[1], uL0[2], uE0[0], uE0[1], uE0[2], uS0[0], uS0[1], uS0[2]] # #    # #   t_begin = 0 #   t_end = 30.7 * Td / T; #      N_plots = 1000 #     step = (t_end - t_begin) / N_plots import scipy.integrate as spi solver = spi.ode(f) solver.set_integrator('vode', nsteps=50000, method='bdf', max_step=1e-6, rtol=1e-12) solver.set_initial_value(y0, t_begin) ts = [] ys = [] i = 0 while solver.successful() and solver.t <= t_end: solver.integrate(solver.t + step) ts.append(solver.t) ys.append(solver.y) print(ts[i], ys[i]) i = i + 1

让我们看看我们得到了什么。从我们选择的起点开始的前29天中，月球的空间轨迹

以及它在黄道平面上的投影。

“嘿叔叔，你给我们什么？！这是圈子！”

首先，它不是一个圆-轨迹的投影从原点到右，下的变化是明显的。其次-什么都没注意到？不，真的吗？

我保证为事实提供理由（基于对帐户错误和NASA数据的分析），即获得的路径偏移不是集成错误的结果。到目前为止，我建议读者相信我的意思-这种位移是月球轨迹受到太阳干扰的结果。再转一圈

什么时候此外，请注意以下事实：根据问题的初始数据，太阳正好位于月亮在每转一圈中移动的一侧。是的，这个无礼的太阳从我们这里窃取了我们挚爱的同伴！哦，是太阳！

我们可以得出结论，太阳引力对月球的轨道影响很大-老妇不会以相同的方式两次向天空行走。半年的移动图片可以（至少在质量上）使人信服（图片可点击）

有意思吗当然可以。天文学通常是一门有趣的科学。

后记

在我学习和工作了近七年的大学-新切尔卡斯克理工学院-每年举行一次北高加索大学理论力学专业的区域奥林匹克竞赛。我们三度参加了全俄奥林匹克运动会。在开幕式上，我们的主要“奥林匹亚”教授A. Kondratenko总是说：“学者Krylov将力学称为精确科学的诗歌。”

我喜欢机械师。多亏了这门科学和我出色的老师，我在生活和事业中所取得的一切成就都实现了。我尊重机械师。

因此，即使他至少是科学博士的三倍和语言学家的四倍，并且已经开发了至少一百万个学习程序，我绝不允许任何人嘲笑该科学并为自己的目的大胆地利用它。我真诚地相信，在流行的公共资源上写文章应该为他们提供仔细的校对，正常的设计（LaTeX公式不是该资源的开发者的心血来潮！）并且没有错误导致导致违反自然法则的结果。后者通常是必须具备的。

我经常告诉我的学生：“计算机可以解放双手，但这并不意味着您也需要关闭大脑。”

为了欣赏和尊重机制，我敦促您，亲爱的读者。我将很乐意回答任何问题，并且如所承诺的那样，我发布了一个示例示例的源代码，该示例在Python中解决了三体问题，我在个人资料中发布了Github 。

感谢您的关注！

动态系统建模：月球如何运动？

引言

1.引力定律以及我们可以从中得出什么结论

2.卫星还是行星？引力范围

3.经典公式中的三体问题

4.方程与无量纲变量的比例分配

5.初始条件的计算：我们获得了初始数据

6.38只鹦鹉和一只鹦鹉翼

7.运动方程式的积分和结果分析

后记

More articles:

动态系统建模：月球如何运动？

引言

1.引力定律以及我们可以从中得出什么结论

2.卫星还是行星？ 引力范围

3.经典公式中的三体问题

4.方程与无量纲变量的比例分配

5.初始条件的计算：我们获得了初始数据

6.38只鹦鹉和一只鹦鹉翼

7.运动方程式的积分和结果分析

后记

More articles:

2.卫星还是行星？引力范围