混沌优化算法求解Kepler方程

2012-10-10杨玲玲

上海理工大学学报 2012年6期

关键词：复杂度次数方程

杨玲玲，马良

（上海理工大学管理学院，上海 200093）

1 问题的提出

1619年Kepler首次推导出了Kepler方程，该方程在物理、数学领域，特别是在经典天体力学研究中都起着重要作用.

在天体力学和轨道计算中，求解Kepler方程

是经常遇到的问题，即在给定平近点角M（0≤M≤π）和偏心率e后去求偏近点角x［1］.在椭圆轨道情况下，0≤e≤1.由于方程中含有超越函数sinx，因此，Kepler方程是个超越方程，这意味着x无法用代数方法求出，只能用数值分析或级数求解，这就要求Kepler方程在实际求解中应尽可能地接近真值.

超越方程常用迭代方法求解，目前应用较多的数值迭代方法是不动点迭代法或牛顿迭代法，但前者收敛速度慢，后者在选取初值上要求较高，而且每次都要计算导数值［2］.

在实际应用中，由于对上千万年的天体轨道作积分需要反复求解Kepler方程，计算量大且耗时，因此，研究快速求解Kepler方程的方法具有重要的意义［3］.本文试图使用一种智能优化算法——混沌优化算法求解Kepler方程，探索从新的角度求解此类问题.

2 现有解法

尽管Kepler方程的提出已近4个世纪，但因其重要性而一直是研究的重要课题.每10～20年就会提出一个新的求解方法.Ng［4］在1979年提出一立方收敛的方法.

1983年Danby和Burkardt提出了具有立方收敛性的Halley方法.

同时，他们基于该方法提出了四次方和五次方的收敛性方法［5］.Colwell，Battin和 Vallado分别提出的方法也是以迭代或级数展开为基础的经典方法［6］.

3 算法设计

混沌运动具有遍历性、随机性及规律性等特点，它能在一定范围内，按其自身规律不重复地遍历所有状态［7］.李兵和蒋慰孙［8］选用Logistic映射做算法

式中，μ为控制参量，取μ＝4.

设0≤a0≤1，n＝0，1，2，…，当μ＝4时，系统完全处于混沌状态，运用载波的方式将混沌状态引入至优化变量中，同时将遍历范围放大到优化变量的取值范围，然后利用混沌变量进行搜索［8］.文献［8］中介绍了算法的基本步骤.

针对Kepler方程设计下面的算法.

Step 1 等价转换式（1）.

将求解式（1）的问题转化为求x，使得目标函数F（x）＝（x－esinx－M）2的值最小.

Step 2 确定x的取值范围.

已知0≤x≤2π，进一步来说，在文献［9］中得出结论

因为f′（x）＝1－ecosx≥0恒成立，0≤e≤1，函数f（x）单调递增，同时

故可得出M≤x≤M＋e.f（x）的图像如图1所示.

图1 f（x）的图形Fig.1 Graph of f（x）

Step 3 初始化.

设置k＝1，a0，最大迭代次数n.

Step 4 根据式（4），产生混沌变量ai.通过式（6），利用载波方式将混沌变量放大到相应优化变量在Step 2中所得到的取值范围.

Step 5 如果F（xi）＜F＊，那么F＊＝F（xi），x＊＝xi；否则，放弃xi.

Step 6 如经若干次搜索后保持不变，则按式（7）进行第二次载波.

式中，x＊为当前最优解；α为较小的调节常数，本文中取小于1的常数.

Step 7 如果F［x（k）］＜F＊，那么F＊＝F［x（k）］，x＊＝x（k）；否则，放弃x（k）.

Step 8 当达到最大迭代次数时，终止搜索，输出最优解x＊，F＊.

不难估算，该算法的时间复杂度为O（n），n为最大迭代次数.

4 数值测试

为检验求解效果，进行了一系列数值测试计算，并与已有方法进行比较.实验所用的硬件为Pentium（R）4CPU 2.80GHz，256MB内存的微机，软件为 Windows XP操作系统和Matlab7.0运行环境.

在现有文献中，大多数计算实例是在M＝0和e＝1附近进行的，因为，在这一点有f′＝f″＝0，在运用传统方法时容易导致解法不收敛，本文首先求解了这类情况，现列出部分结果.表1和表2分别列出了迭代次数为1 000和10 000时，在临界情况（零平近点角和等于1的偏心率，即M＝0，e＝1）附近的测试结果，混沌初值a0＝0.201 7.