整整半个月过去了，怎么一直无人对此有兴趣？
说白了，我就是想用数值积分方法精确计算水星近日点进动值变化规律，从而进一步证明百年前西蒙·纽康的计算值错误。
纽康当年的计算值，是根据1855年Le Verrier 发表的7阶Lagrange方程改进展开式计算得到的。这改进展开式其实仍有不少错误，其中一个大错误直到1985年才被Murray纠正，可见当年纽康计算的牛顿引力每百年摄动531.9角秒值根本就不值得可信，精确值究竟是多少？我想非用Mathematica编程进行数值积分不可。

许剑伟设计的“寿星天文历”软件的确非常不错，但是他介绍《Astronomical Ephemeris》等权威天文历书中提出的“力学时”与天体力学概念似乎有点混淆，他们把地球椭圆轨道速度的变化误解为(广义)“相对论效应”，并引入他的寿星软件中被牧夫天文论坛的高手们佩服得五体投地。
其实他所译著的《天文算法》非常简单，用的是半解析经验公式计算，根本就没有用上科威耳积分公式进得数值积分计算，与牛顿平方反比力作用过程也没有直接关系，所以我怀疑历史上(法国)由《天文算法》误差得到的行星轨道近日点进动等计算数据并不可靠。

决定重新对Lagrange摄动方程尝试30级的级数展开计算。当年西蒙·纽康用手工计算的Le Verrier 的7阶展开式错误是难免的，而且那种平近点角展开方法计算量太大精度太低(至今用电脑程序也还只得到一个8阶式)。以金星与水星摄动计算为例，7阶只能达到两位有效值精度，而30阶可达10位有效摄动值。
30阶级数展开只展开到真近点角的指数级数，因为在Mathematica中可以用反函数序程包把真近角直接计算到任意所需精度。
完成30阶级数展开(扁心率展开到12级)，预期符号运算总项数在10万以内，展开后直接进行符号运算普通电脑的内存可能太小，可以分级打包输入数据实现数值计算。例如把角度向量化为一个四项式，直接对这四项式的30次方用电脑进行符号积分运算几乎是不可能的，转成数值系数计算就简单了。
先完成30阶的全部级数展开计算，输入程序包保存，再取出输入全部轨道根数，激活就能得计算结果，还可以找出一些新的长周期共振摄动项。
模拟行星真实运动状态的数值积分问题，因为我不懂自动迭代运算编程，只好以后有机会再说。

走了很多弯路，先是把勒让德级数展开次序弄错得到不规则的发散结果。勒让德级数展开到30级有256项，除前几十项简单点外，其它所有项都可以再展开由六重数列组合而成的复杂系统，至少有上百万个数据需要计算。之前分项输入工作量太大太大，而且容易出错，现在改用函数输入就简单了，但计算机的负荷太大，几百万数据一次计算需二十多分钟，完成全部计算量要上百次这种计算，计算机要改朝换代了，过几天再买台高性能电脑回。
计算行星轨道牛顿引力摄动，用拉格朗日摄动公式计算初步得出的水星轨近日点进动值根本不是那么回事了，每百年竟高达1529.493"，而两年前我用Excel 表格计算的结果则是528.767" (全部逐项展开估丢了很多项数)，这次用Mathematica 软件计算的应该不会，程序设计和输也应该没错，究竟错在哪还需进一步。
对几百万全都展的数据计算平均值，有很多很多技巧，一般人学不过来，其中关键的是阶乘和双阶乘的巧妙应用，而且为了保证计算精度展开到三十级是完全必要，百多年前天文学家用手工计算只展开到7级，在现代计算机应用手段面前实在太小儿科了。

水星轨道的牛顿引力摄动近日点进动值，完整的计算结果应该是每百年588"左右，而不是百多年前天文学家们计算的532"。现在的观测值是573.57"，超出的14"多则是我的理论修正值。就是因为这14"多才确保了引力系统的物质守恒和黑洞不存在以及宇宙熵平衡的前提条件，从而促使引力质量不断向惯性质量转化而温度升高辐射加强，最后引力物质转换成光的辐射能全部被带走而弥散到整个宇宙空间。
注意：光波是具有完整的物质质量属性的，是物质存在的第二种形式，第三种形式则是充满整宇宙空间以太和引力场。

换了AMD 8核CPU和技嘉970主版的运算速度居然比之前4核的还慢了！
经多台电脑反复测试，同一数据计算时间，双核、4核和8核的分别是25分钟、18分钟和21分钟，CPU的运行使用率分别是50%、25%和13%，可能Wolfram Mathematica 序程是按单核CPU满负荷设计的，所以增加CPU内核数不可能加快运算速度，但与硬盘质量和其它配置密切相关(我的新电脑用的是金士顿固态硬盘)。
过几天再逐步公布计算结果。

八大行星的牛顿引力摄动近日点和升交点进动值，用Lagrange摄动方程理论计算结果分别(括号内是观测值，单位为每百年角秒)是：
近日点进动值：水星 528.2638(573.57)，金星 346.2876(-108.80)，地球 1456.9983(1198.28)，火星 1552.0332(1560.78)，木星 643.7413(839.93)，土星 2144.8830(-1948.89)，天王星 292.3872(1312.56)，海王星 65.8545(-844.43)。
升交点进动值：水星 -528.1635(-446.30)，金星 -994.2910(-996.89)，地球 -5628941.6785(-18228.25)，火星 -1084.4419(-1020.19)，木星 671.2468(1217.17)，土星 -958.8588(-1591.05)，天王星 276.5821(-1681.40)，海王星 -22.0242(-151.25)。

最后计算结果还是与我五年前用 Excel 表格计算结果非常接近，无论是高斯公式还是拉格朗日公式多次计算结果都一样，用Mathematica 的精确计算是(每百年角秒)：
近日点进动值：水星 528.7694，金星 12.0079，地球 1135.3418，火星 1593.9554，土星 643.0834，木星 1658.2315，天王星 311.0234，海王星 71.1640。
升交点进动值：水星 -451.4844，金星 -994.2752，地球 -5754713.4，火星 -1059.5182，土星 645.1411，木星 -928.1120，天王星 279.8878，海王星 -23.1862。

大太神奇了，把一个摄动方程设置几百个关于时间的不同反函数和隐函数构成的数值积分，居然可以直接计算并轻易完成高精度计算量，之前我简直不敢相信！现在不知能否进一步直接计算这类微分方程。
遗憾的是用这种方法计算轨道摄动虽可以反映某一阶段真实状态，却很难得到平均值，也许再用作图方法可以找出平均波动趋势。

(上楼图片发错了，删除重发这一张)这是从历元J2000开始100年内水星轨道近日点数据分布状况，这种波动是轨道半长轴在土星方向的表现，具有典型的约11年周期性，单位是每天角秒。

用Mathematica8计算原数据中因海王星的轨扁心率输入有误，行星轨道牛顿引力摄动重新计算后所得最后结果为：
近日点每百年进动值(括符内为观测值)：水星528.769"(573.57")，金星10.934"(-108.80")，地球1143.22(1198.28")，火星1593.54"(1560.78")，土星641.59"(839.93")，木星1660.66"( -1948.89")，天王星291.435"(1312.56")，海王星82.626(-844.43")。
升交点每百年进动值(括符内为观测值)：水星-451.484"(-446.30")，金星-1000.55"( -996.89")，地球-5792934.37"(-18228.25")，火星-1062.37"(-1020.19")，土星645.147"(1212.17")，木星-939.924"(-1591.05")，天王星270.448"(-168140")，海王星 -21.530(-151.25")。
因金星与地球轨道半径比值较大，用拉格朗日摄动公式展开到30级的计算精度还不够，同时有7重级数展开计算量太大电脑速度太慢，故采用高斯摄动公式20次迭代运算，所以重新计算轻松得到了足够高精度的计算结果。
相关人员不应再沉迷于水星问题的43"了，看看其它行星问题又如何？不过这里还有太阳扁率摄动和银河系自转数据未加入修正，是否还有其它更大因素对所有行星的影响不得而知。