楊玲玲， 馬 良

（上海理工大學(xué) 管理學(xué)院，上海 200093）

1 問題的提出

1619年Kepler首次推導(dǎo)出了Kepler方程，該方程在物理、數(shù)學(xué)領(lǐng)域，特別是在經(jīng)典天體力學(xué)研究中都起著重要作用.

在天體力學(xué)和軌道計算中，求解Kepler方程

是經(jīng)常遇到的問題，即在給定平近點角M（0≤M≤π）和偏心率e后去求偏近點角x［1］.在橢圓軌道情況下，0≤e≤1.由于方程中含有超越函數(shù)sinx，因此，Kepler方程是個超越方程，這意味著x無法用代數(shù)方法求出，只能用數(shù)值分析或級數(shù)求解，這就要求Kepler方程在實際求解中應(yīng)盡可能地接近真值.

超越方程常用迭代方法求解，目前應(yīng)用較多的數(shù)值迭代方法是不動點迭代法或牛頓迭代法，但前者收斂速度慢，后者在選取初值上要求較高，而且每次都要計算導(dǎo)數(shù)值［2］.

在實際應(yīng)用中，由于對上千萬年的天體軌道作積分需要反復(fù)求解Kepler方程，計算量大且耗時，因此，研究快速求解Kepler方程的方法具有重要的意義［3］.本文試圖使用一種智能優(yōu)化算法——混沌優(yōu)化算法求解Kepler方程，探索從新的角度求解此類問題.

2 現(xiàn)有解法

盡管Kepler方程的提出已近4個世紀，但因其重要性而一直是研究的重要課題.每10～20年就會提出一個新的求解方法.Ng［4］在1979年提出一立方收斂的方法.

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

同時，他們基于該方法提出了四次方和五次方的收斂性方法［5］.Colwell，Battin和 Vallado分別提出的方法也是以迭代或級數(shù)展開為基礎(chǔ)的經(jīng)典方法［6］.

3 算法設(shè)計

混沌運動具有遍歷性、隨機性及規(guī)律性等特點，它能在一定范圍內(nèi)，按其自身規(guī)律不重復(fù)地遍歷所有狀態(tài)［7］.李兵和蔣慰孫［8］選用Logistic映射做算法

式中，μ為控制參量，取μ＝4.

設(shè)0≤a0≤1，n＝0，1，2，…，當μ＝4時，系統(tǒng)完全處于混沌狀態(tài)，運用載波的方式將混沌狀態(tài)引入至優(yōu)化變量中，同時將遍歷范圍放大到優(yōu)化變量的取值范圍，然后利用混沌變量進行搜索［8］.文獻［8］中介紹了算法的基本步驟.

針對Kepler方程設(shè)計下面的算法.

Step 1 等價轉(zhuǎn)換式（1）.

將求解式（1）的問題轉(zhuǎn)化為求x，使得目標函數(shù)F（x）＝（x－esinx－M）2的值最小.

Step 2 確定x的取值范圍.

已知0≤x≤2π，進一步來說，在文獻［9］中得出結(jié)論

因為f′（x）＝1－ecosx≥0恒成立，0≤e≤1，函數(shù)f（x）單調(diào)遞增，同時

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

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

Step 3 初始化.

設(shè)置k＝1，a0，最大迭代次數(shù)n.

Step 4 根據(jù)式（4），產(chǎn)生混沌變量ai.通過式（6），利用載波方式將混沌變量放大到相應(yīng)優(yōu)化變量在Step 2中所得到的取值范圍.

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

Step 6 如經(jīng)若干次搜索后保持不變，則按式（7）進行第二次載波.

式中，x＊為當前最優(yōu)解；α為較小的調(diào)節(jié)常數(shù)，本文中取小于1的常數(shù).

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

Step 8 當達到最大迭代次數(shù)時，終止搜索，輸出最優(yōu)解x＊，F(xiàn)＊.

不難估算，該算法的時間復(fù)雜度為O（n），n為最大迭代次數(shù).

4 數(shù)值測試

為檢驗求解效果，進行了一系列數(shù)值測試計算，并與已有方法進行比較.實驗所用的硬件為Pentium（R）4CPU 2.80GHz，256MB內(nèi)存的微機，軟件為 Windows XP操作系統(tǒng)和Matlab7.0運行環(huán)境.

在現(xiàn)有文獻中，大多數(shù)計算實例是在M＝0和e＝1附近進行的，因為，在這一點有f′＝f″＝0，在運用傳統(tǒng)方法時容易導(dǎo)致解法不收斂，本文首先求解了這類情況，現(xiàn)列出部分結(jié)果.表1和表2分別列出了迭代次數(shù)為1 000和10 000時，在臨界情況（零平近點角和等于1的偏心率，即M＝0，e＝1）附近的測試結(jié)果，混沌初值a0＝0.201 7.

表1 最大迭代次數(shù)為1 000時的結(jié)果Tab.1 Results of n＝1000

從表1和表2中可以看出，迭代次數(shù)越大，求解精度越好.由于算法的時間復(fù)雜度較低，算法的運行時間基本可以忽略不計.

文獻［10］中指出，用改進的牛頓法求解時可能收斂至一個或多個錯誤的值，例如，當M＝1.347，e＝0.66，迭 代 序 列 在 0.448 902，0.354 721，0.469 946和0.362 751這4個值上循環(huán)，而此時的正確解是1.958 11［10］.所以，本文也對該情況進行了求解，求解結(jié)果為

表2 最大迭代次數(shù)為10 000時的結(jié)果Tab.2 Results of n＝10 000

本文還測試了文獻［3］和文獻［10］中的一些算例，結(jié)果如表3所示.

表3 部分算例測試結(jié)果Tab.3 Part of test results

測試結(jié)果表明，混沌優(yōu)化算法對于Kepler方程求解的有效性，該算法較低的時間復(fù)雜度提升了運算的效率，且免除了多次求導(dǎo)的復(fù)雜運算，避免了在臨界狀態(tài)容易不收斂的缺點.

5 結(jié) 論

常見的求解Kepler方程的方法是簡單迭代法、牛頓法等經(jīng)典方法，本文運用智能優(yōu)化方法——混沌優(yōu)化算法求解該問題.算法的時間復(fù)雜度較低，具有較高的運算效率，同時避免了傳統(tǒng)方法中多次求導(dǎo)等復(fù)雜計算，且取得了滿意的結(jié)果.對于Kepler方程，本文所給出的混沌優(yōu)化算法是一種簡單而又快速地獲得有效解的途徑.

［1］周慶林，黃天衣.Kepler方程的解法［J］.天文學(xué)報，1988，29（1）：106－112.

［2］束雄英，李紅.Kepler方程的一種迭代加速［J］.航空計算技術(shù)，2005，35（1）：41－44.

［3］王玉詔，鐘雙英，孫威，等.Kepler方程的六階迭代解法［J］.江西科學(xué)，2009，27（6）：790－792.

［4］Colwell P.Solving Kepler’s equation over three centuries［M］.Richmond：Willmann－Bell，1993.

［5］Danby J M A，Burkardt T M.The solution of Kepler’s equation［J］.Celestial Mechanics，1983，31（2）：317－328.

［6］Davis J J.Sequential solution to Kepler’equation［J］.Celest Mech Dyn Astr，2010，108（7）：59－72.

［7］苗東升.系統(tǒng)科學(xué)大學(xué)講稿［M］.北京：中國人民大學(xué)出版社，2007.

［8］李兵，蔣慰孫.混沌優(yōu)化方法及其應(yīng)用［J］.控制理論與應(yīng)用，1997，14（4）：613－615.

［9］Smith G R.A simple，efficient starting value for the iterative solution of Kepler’s equation［J］.Celestial Mechanics，1979，19（2）：163－166.

［10］岳錦海，黃天衣.橢圓Kepler方程求解中的非線性現(xiàn)象［J］.南京大學(xué)學(xué)報（自然科學(xué)版），1998，34（1）：21－28.