物理科学探疑-网友天空-宇宙的观念

月球轨道计算

芦光耀

    计算机可精确地求解星球运动微分方程。本算法适用于用计算机求解任意多体问题的星球运动微分方程，也能精确地确认万有引力公式本身和实际情况的相差程度。

一．计算机计算星球轨道原理。     以地球绕日运行为例，如图1-1。地球受力如图1-2。     为叙述方便，以下用矢量形式表达。

    因为其中M_太为太阳质量，M_地为地球质量，G为万有引力恒量，a为地球加速度,为单位矢量。所以我们可按等时间间隔dt（即等步长），以微分形式从地球的初值点逐点向下推算。设t=0时，地球的初值点为r₀,v₀和于是，地球经dt时间从初值点到达第一点,递推式为,

    由于dt是人为设定的，是已知的，因此地球到达1点的近似值v,r和a可由上式算出，算出1点值后，可把1点值作为初值，按步长dt继续推算出下一点的值，如此，可推算到第n点。由于dt值取得越小，递推的精度越高，我们可据此来控制计算误差。

    设要计算地球在t=T时的r值，要求计算误差为e，t=0时的初值r₀,a₀,v₀为已知。我们可将0到T的时间间隔划分为n个dt，即令计算步长dt=T/n,然后根据上述，按步长dt从t=0时的初值点推算到t=n·dt=T时的r值。然后将dt二分，即令计算步长dt₁=dt/2,再按此新步长值dt₁从t=0时的初值点算到t=2·n·dt₁=T时的r值为r₂，比较一下二分前后的r值，即看一看是否满足条件r₂ - r<e,若满足条件，则r₂的值即为所要求的r值，若不满足条件，则继续二分，按新的步长值再从t=0时的初值点算到t=T时刻，直到新算出的r值满足条件r₂ - r<e，其中r和r₂分别为二分前后的r值。注意，上面所说的矢量比较包含其大小比较和其方向比较。

    以上为矢量表达，实际计算中，可将矢量v,r，a分别在x,_y轴上投影，可得v_x,v_y,r_x,r_y,a_x,a_y。于是，它们的初始值分别为v_0x,v_0y,r_0x,r_0y,a_0x,a_0y。下面给出地球从初值点经dt时间运行到下一点的递推式，

    控制误差范围的条件为分别为二分前后在x,y方向的r值。

    以上仅为计算机计算地球轨道的原理。实际上，每二分一次，从0到T时间范围内的dt数量将增加一倍，计算机计算的工作量也将增加一倍。由于计算机的计算速度有限，因此二分次数也是有限的。为提高计算精度，减少计算机的计算工作量，有一些标准化的方法（注1），在此不再熬述。

    由上述可知，计算机计算星球轨道主要有两个要点。一是列出递推式，二是确定误差范围的条件。

    月球轨道计算见下页。

注释1：参见“计算机数值计算方法及程序设计”一书。该r书由周煦编著。于2004年10月由机械工业出版社出版。

二．月球轨道计算

    由于地球的运动直接影响月球的运动，因此，先来分析一下地球的受力，如图1-3所示。

    在图1-3中，o₂x₂y₂z₂坐标系是动坐标系，原点在地球中心。该坐标系跟随地球作平动，且三个坐标轴x₂,y₂,z₂始终分别平行于x,y,z三个坐标轴。r₁ 是地球的位置矢量，r是月球的位置矢量, r₂ 是月球相对地球的位置矢量。

    F月地是月球对地球的引力，F太地是太阳对地球的引力。设r₁ 与x,y,z轴的夹角分别为α₁，β₁，γ₁，r与x,y,z轴的夹角分别为α，β，γ，r₂ 与x₂,y₂,z₂轴的夹角分别为α₂,β₂,γ₂,则，地球在x,y,z方向所受合力为：

    因此，地球在x,y,z方向的加速度：

    月球的受力如图1-4所示。月球在x,y,z方向所受合力为:

    其中，F_太月为太阳对月球的引力，F_地月为地球对月球的引力。因此，月球的加速度为：

设a的初值为的初值为这样，地球和月球从各自的初值点同时出发，经dt时间后，地球就到达了它的下一点于是可得如下递推式：

(见下页)

控制计算误差的6个条件为：

    其中分别为二分前后算出的地球坐标。再次说明一下，以上月球轨道的计算仅是计算机计算原理，实际编程应采取一些标准化方法，以提高计算精度，减少计算机的计算工作量。

    目前，在月球轨道计算上，我已做到了，一天的计算误差e<0.001米（即在x,y,z轴方向的计算误差e）,也就是说一年的计算误差e<365×0.001=0.365米。要核实万有引力公式本身和实际情况的相差程度，可取两组实际观测值，一组观测值作为计算的初值，另一组观测值作核实之用，即核实用万有引力公式来计算的星球轨道的准确程度。下面采用一组实际观测值（注2）作为计算初值，让计算机来计算一下月球的轨道。初值为：

    下面给出计算结果。

    以上计算的时间范围是2006年。月球轨道半径最大和最小值都是指平均值。观测值由紫金山天文台的工作人员提供。以上计算取计算时间t=366天，坐标计算误差e< 0.366米。计算周期时，时间计算误差小于0.05秒。由于天文台提供的月球数据是相对地球坐标系的，地球坐标系和本文所述的动坐标系的关系是，将地球坐标系的yz平面以x轴为转轴旋转一个黄赤交角，就是本文所述的动坐标系。本文取黄赤交角为23^o26′。

注释2：作为计算初值的观测值的对应时刻为，2006年北京时间3月15日7时47.5分。该时刻恰好为半影月食的食甚时刻。为方便计算， 本文把该时刻定为零时刻。

下面给出从2006年3月月食食甚时刻计算到9月月食食甚时刻的地球和月球坐标。数据如下。
地球坐标（单位：米）

作者：芦光耀
日期：2007年5月4日。
住址：北京市丰台区东高地三角地37栋丁门四楼111号。
邮编：100076
电话：68751647
E-mail：hulahuaer@hotmail.com

     版权所有，保留一切权力，未经授权使用将追究法律责任 版权说明  © Copyright  Authors
物理科学探疑

返回首页