吳連大 張明江

(1 中國科學(xué)院紫金山天文臺南京 210023)

(2 中國科學(xué)院空間目標(biāo)與碎片觀測重點實驗室南京 210023)

1 引言

在天體力學(xué)的分析理論中,需要將攝動函數(shù)展開為時間和軌道根數(shù)的顯函數(shù),這涉及到兩個最常用的特殊函數(shù):傾角函數(shù)和Hansen系數(shù)[1].在衛(wèi)星動力學(xué)中,地球引力場和日月第3體引力攝動的攝動函數(shù)展開即與Hansen系數(shù)有關(guān)[2–3].

Balmino[4]、Wnuk[5]、Giacaglia[6]、McClain[7]和Proulx等[8]給出了Hansen系數(shù)的各種直接計算表達(dá)式.利用這些Hansen系數(shù)的直接計算表達(dá)式,對軌道偏心率e直接求導(dǎo),我們可以得到相應(yīng)的Hansen系數(shù)導(dǎo)數(shù)的直接計算表達(dá)式.基于Hansen系數(shù)的定積分表達(dá)式及其直接求導(dǎo)的導(dǎo)數(shù)表達(dá)式,也可以實現(xiàn)Hansen系數(shù)及其導(dǎo)數(shù)的直接計算.此外,Wnuk[5]還給出了另一種Hansen系數(shù)導(dǎo)數(shù)的直接計算表達(dá)式.這樣,對于Hansen系數(shù)及其導(dǎo)數(shù)的直接計算,本文回顧總結(jié)出7種主要方法:Balmino方法[4]、Wnuk方法[5]、Wnuk直接求導(dǎo)方法、Giacaglia方法[6]、兩種McClain方法[7–8]和定積分方法.這些計算方法涉及的具體數(shù)學(xué)表達(dá)式,參見附錄.需要說明的是,附錄相關(guān)計算方法的數(shù)學(xué)表達(dá)式中,求和指標(biāo)的符號遵照了原文獻(xiàn)的形式,以便讀者查閱原文.

本文針對這7種Hansen系數(shù)及其導(dǎo)數(shù)的直接計算方法,重點探討各種方法的計算效率以及穩(wěn)定性,進(jìn)而得到一些有益的結(jié)論,為實際應(yīng)用中Hansen系數(shù)及其導(dǎo)數(shù)的計算提供參考.

2 計算效率

為了比較分析,我們編制了這7種Hansen系數(shù)及其導(dǎo)數(shù)的直接計算方法的Fortran程序,包括:

(1) Balmino方法[4],即附錄中表達(dá)式(1a)和(1b),其中指標(biāo)t的求和限為0–50;

(2) Wnuk方法[5],即附錄中表達(dá)式(2a)和(2b),其中指標(biāo)s的求和限為0–70,指標(biāo)t的求和限為?70至70;

(3) Wnuk直接求導(dǎo)方法,即附錄中表達(dá)式(3a)和(3b),其中指標(biāo)s的求和限為0–70,指標(biāo)t的求和限為?70至70,注意附錄中表達(dá)式(2a)和(3a)實質(zhì)上是相同的;

(4) Giacaglia方法[6],即附錄中表達(dá)式(4a)和(4b),其中,當(dāng)指標(biāo)時,指標(biāo)s的求和限為0–70;當(dāng)指標(biāo)0時,指標(biāo)t的求和限為0–70;

(5) McClain方法1[7–8],即附錄中表達(dá)式(5.1a)和(5.1b),其中指標(biāo)i的求和限為0–70;

(6) McClain方法2[7–8],即附錄中表達(dá)式(5.2a)和(5.2b),其中指標(biāo)i的求和限為0–70;

(7)定積分方法,即附錄中表達(dá)式(6a)和(6b),使用分段兩點Gauss積分方法,分段數(shù)為2000.

利用上述7種方法,計算偏心率函數(shù)Glpq=(l∈[2,30],p∈[0,l],q∈[?2,2],e=0.1)及其導(dǎo)數(shù),各種方法的計算時間,如表1所示.計算設(shè)備的基本配置參數(shù)為:處理器Intel?Xeon?W-2125 CPU @ 4.00 GHz (8 CPUs),主頻約4.0 GHz;內(nèi)存32768 MB RAM.由表1可見,各種方法的計算效率差別很大,在具體使用時需要注意.另外,需要注意的是,目前本文中各種計算方法所采用的相應(yīng)指標(biāo)的求和限未必最優(yōu),讀者在使用時應(yīng)根據(jù)需要進(jìn)行必要的數(shù)值試驗,以確定指標(biāo)合適的求和限.實際上,各種計算方法中指標(biāo)的求和限選取,主要取決于偏心率的大小和需要計算的Hansen系數(shù)的階次.

表1 各種方法的計算時間Table 1 Computing time of various methods

3 穩(wěn)定性分析

Wnuk[5]對相關(guān)Hansen系數(shù)及其導(dǎo)數(shù)直接計算方法的穩(wěn)定性,有過簡約的評價:

(1)利用改造Bessel函數(shù)(漸近展開)算法的Wnuk方法,對應(yīng)附錄(2a)和(2b)式,是數(shù)字穩(wěn)定的,對于偏心率e<1以及高階l、p和q,能進(jìn)行偏心率函數(shù)Glpq=及其導(dǎo)數(shù)的快速計算;

(2)Kaula方法[2](與McClain方法2本質(zhì)上相同,對應(yīng)附錄(5.2*)、(5.2a)和(5.2b)式),對于大偏心率是數(shù)字不穩(wěn)定的;

(3)定積分方法,對應(yīng)附錄(6a)和(6b)式,對于大偏心率e情形計算得很好,但由于計算時間和數(shù)字不穩(wěn)定性,在實際應(yīng)用中對高階l、p和q的偏心率函數(shù)計算受到限制.

下面,我們結(jié)合上述7種Hansen系數(shù)及其導(dǎo)數(shù)的直接計算方法的Fortran程序計算算例,較深入地分析這些方法的計算穩(wěn)定性.

3.1 Hansen系數(shù)計算的不穩(wěn)定性例子

計算結(jié)果表明:對于小偏心率(e <0.2),各種方法計算均是穩(wěn)定的.但是,對于大偏心率,Hansen系數(shù)計算會出現(xiàn)不穩(wěn)定的情況.例如,對于偏心率e=0.6,Hansen系數(shù)的相關(guān)計算結(jié)果,見表2.

由表2中計算結(jié)果可見,對于Hansen系數(shù)的直接計算,Wnuk方法(Wnuk直接求導(dǎo)方法)和Giacaglia方法(即附錄(3a)和(4a)式)的計算結(jié)果符合得較好(附錄(2a)與(3a)式實質(zhì)相同);定積分方法(即附錄(6a)式)的計算結(jié)果與Wnuk方法、Giacaglia方法整體上一致.然而,Balmino方法和McClain方法1 (即附錄中(1a)和(5.1a)式)的計算結(jié)果與Wnuk方法相比,從q≥?4開始就符合得不好了;McClain方法2 (即附錄中(5.2a)式)的計算結(jié)果與Wnuk方法相比,從q≥6開始符合得不好了.計算表明:這些方法在大偏心率高階的情況下,出現(xiàn)了不穩(wěn)定的情況.這種不穩(wěn)定現(xiàn)象的根源在于計算機(jī)字長不夠,如果采用4精度計算,結(jié)果就穩(wěn)定了.

表2 Hansen系數(shù)計算算例(e=0.6)Table 2 The calculation example of Hansen coefficients (e=0.6)

表2 續(xù)Table 2 Continued

3.2 計算方法穩(wěn)定性的判別

為了比較各種方法的優(yōu)劣,需要一個判別計算方法穩(wěn)定性的準(zhǔn)則.實際上,利用直接計算方法計算得到的Hansen系數(shù),通過判別這些數(shù)據(jù)是否滿足Hansen系數(shù)的遞推關(guān)系,就可判別方法的穩(wěn)定性.

不難驗證,對于小偏心率,各種方法的計算結(jié)果,均滿足遞推關(guān)系.但是,對于大偏心率,就不一定滿足了.舉一個例子,偏心率函數(shù)的遞推關(guān)系式為[6]:

對于偏心率e=0.75,指標(biāo)l=30、p=29和q=?1,Wnuk方法和McClain方法1,即附錄(2a)和(5.1a)式程序計算結(jié)果,見表3.要判別Wnuk方法和McClain方法1的計算結(jié)果是否穩(wěn)定,只需將計算結(jié)果代入上述遞推關(guān)系式即可.對于Wnuk方法,此時Glpq=0.2056355791128×103、Gl?1,p,q+1=0.4210025874464×102、Gl?1,p?1,q?1=0.9079971968072×102、Gl?2,p?1,q=0.1542231019741×102;而l=30、l?2p=?28、l?2p+q=?29、l?1=29,于是遞推關(guān)系式右邊=94.5116409715.42231020 (遞推關(guān)系式左邊),顯然這個結(jié)果不能認(rèn)為是穩(wěn)定的.對于McClain方法1,通過類似判別,相應(yīng)的計算結(jié)果更不穩(wěn)定.

表3 穩(wěn)定性判別數(shù)據(jù)(雙精度計算)Table 3 The data for stability discrimination (double precision computation)

進(jìn)一步分析表明,Wnuk方法的計算結(jié)果不穩(wěn)定是由于求和范圍不夠,而McClain方法1不穩(wěn)定是由于計算機(jī)字長不夠.我們擴(kuò)大Wnuk方法的求和范圍,在McClain方法1中利用4精度計算,計算結(jié)果如表4所示.由表4可知,這兩種結(jié)果可以認(rèn)為是穩(wěn)定的.應(yīng)該說明,由于McClain方法1的計算效率很高,即使采用4精度來計算時間也較短,可以滿足實用要求.

表4 穩(wěn)定性判別數(shù)據(jù)Table 4 The data for stability discrimination

3.3 小結(jié)

根據(jù)7種Hansen系數(shù)及其導(dǎo)數(shù)直接計算方法的Fortran程序計算情況,結(jié)合上述分析,可以對Hansen系數(shù)計算,得出如下初步結(jié)論:

(1)對于小偏心率(e<0.2)軌道,各種計算方法均可以滿足計算精度的要求;

(2)對于大偏心率軌道,Wnuk方法(Wnuk直接求導(dǎo)方法)和Giacaglia方法計算結(jié)果較好;

(3)可以利用遞推關(guān)系來判別計算結(jié)果是否穩(wěn)定;

(4)Hansen系數(shù)計算的主要困難是計算機(jī)字長不夠.

Hansen系數(shù)計算不穩(wěn)定的情況,均出現(xiàn)在大偏心率軌道情形,而小偏心率軌道沒有問題.對于大偏心率軌道,分析Hansen系數(shù)計算丟失有效位數(shù)的原因,主要是由于Hansen系數(shù)計算,均是由級數(shù)求和得到,即y=∑xi.由于xi中含有階乘和二項式系數(shù),在偏心率較大時,其數(shù)量級相差很大,求和時會損失計算精度.最糟糕的是在求和時數(shù)量級大的數(shù)正負(fù)相消,余下部分(即為Hansen系數(shù)計算結(jié)果)是一個小數(shù),如果這個小數(shù)比最大數(shù)的末位數(shù)值還小,計算結(jié)果就沒有有效數(shù)字,不穩(wěn)定問題就呈現(xiàn)出來.

由于大多數(shù)人造衛(wèi)星采用小偏心率軌道,使用無奇點根數(shù)的攝動計算是必須的.針對無奇點根數(shù)的攝動計算,我們需要的是Hansen系數(shù)核,即須從Hansen 系數(shù)中提取e|k?m|因子,同時Hansen系數(shù)核本身不存在偏心率e為零的小分母問題.Wnuk方法,即附錄(2a)或(3a)式,不能滿足這一要求,不能計算Hansen系數(shù)核及其導(dǎo)數(shù),因而不適用于無奇點根數(shù)的攝動計算.采用4精度計算的McClain方法1,即附錄(5.1a)式,可以兼顧大小偏心率的Hansen系數(shù)計算,也許是較好的選擇.

4 結(jié)論

本文回顧總結(jié)了7種Hansen系數(shù)及其導(dǎo)數(shù)的直接計算方法,給出了相關(guān)方法的直接求導(dǎo)的導(dǎo)數(shù)表達(dá)式,比較分析了這些方法的計算效率和計算穩(wěn)定性.研究表明:Hansen系數(shù)的遞推關(guān)系可以用來判別計算結(jié)果的穩(wěn)定性.值得指出的是,Wnuk方法(雙精度計算)和McClain方法(4精度計算)是穩(wěn)定的,可以用來計算人造衛(wèi)星軌道攝動.由于

大多數(shù)人造衛(wèi)星采用小偏心率軌道,需要計算無奇點攝動,推薦使用McClain方法1 (4精度計算),即附錄中(5.1a)和(5.1b)式.

附錄

A.1 Balmino方法

基于Balmino給出的Hansen系數(shù)表達(dá)式[4],利用楊輝三角,我們得到下列Hansen系數(shù)的計算表達(dá)式:

上述表達(dá)式(1a)和(1b)中,s必須恒為正.如果s <0,則利用Hansen系數(shù)的對稱性進(jìn)行計算,即表達(dá)式中的m和s均改變符號.

A.2 Wnuk方法

Wnuk給出下列Hansen系數(shù)及其導(dǎo)數(shù)的計算表達(dá)式[5]:

其中,

Jt(ke)是關(guān)于ke的Bessel函數(shù).

Wnuk建議,按如下方法計算Bessel函數(shù)Jt(ke)[5,9]:

值得指出的是:J0(ke)有許多漸近展開式,上式中ξ和η的計算表達(dá)式引自文獻(xiàn)[10],試算表明,這樣計算J0(ke)的相當(dāng)精度約為10?8,已能基本滿足計算Hansen系數(shù)的要求;pt的遞推計算,參見文獻(xiàn)[9],從pt+1=0開始遞推,取t=2N,N為需要計算的Jt(ke)的指標(biāo)t的最大值;當(dāng)ke=0時,上式中pt的計算有0分母,需要特殊處理.

當(dāng)t<0或ke<0時,需用下式計算Bessel函數(shù)[9]:

A.3 Wnuk直接求導(dǎo)方法

Hansen系數(shù)的計算表達(dá)式(2a)[5],可改寫為如下形式:

其中,m前的正負(fù)號可以根據(jù)k?t?m是否大于0來決定,

其中,

A.4 Giacaglia方法

其中,

我們將(4a)式中偏心率函數(shù)Glpq,直接對軌道偏心率e求導(dǎo),得到相應(yīng)的偏心率函數(shù)導(dǎo)數(shù)的計算表達(dá)式:

其中,

A.5 兩種McClain方法

(1) McClain方法1

我們將(5.1*)式改寫為下列形式,并將其對軌道偏心率e直接求導(dǎo),有:

其中,

(2) McClain方法2

我們將(5.2*)式改寫為下列形式,并將其對軌道偏心率e直接求導(dǎo),有:

A.6 定積分方法

利用軌道偏近點角E作為積分變量,Hansen系數(shù)及其導(dǎo)數(shù)的定積分計算表達(dá)式為:

其中,r為目標(biāo)向徑,a、e、M和f分別為軌道半長徑、偏心率、平近點角和真近點角.