鄧子辰　李慶軍

西北工業(yè)大學(xué)工程力學(xué)系, 西安 710072; ? E-mail: dweifan@nwpu.edu.cn

?

精細(xì)指數(shù)積分法在衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)中的應(yīng)用

鄧子辰?李慶軍

西北工業(yè)大學(xué)工程力學(xué)系, 西安 710072; ? E-mail: dweifan@nwpu.edu.cn

編隊(duì)飛行衛(wèi)星間的距離遠(yuǎn)小于衛(wèi)星的軌道半徑, 其動(dòng)力學(xué)方程表現(xiàn)為弱非線(xiàn)性。針對(duì)弱非線(xiàn)性方程的求解, 提出精細(xì)指數(shù)積分方法, 用精細(xì)積分法求解指數(shù)積分方法中的指數(shù)矩陣。用精細(xì)指數(shù)積分法和 Runge-Kutta方法, 在不同條件下求解弱非線(xiàn)性方程的算例, 驗(yàn)證了精細(xì)指數(shù)積分法的有效性。通過(guò)Lagrange方程,建立衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)方程的半線(xiàn)性形式, 用精細(xì)指數(shù)積分方法與 Runge-Kutta 方法求解方程。數(shù)值計(jì)算結(jié)果表明, 與同階的 Runge-Kutta 求解弱非線(xiàn)性微分方程相比, 精細(xì)指數(shù)積分法具有更高的精度, 為衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)仿真提供了一種有效的數(shù)值算法。

指數(shù)積分方法; 精細(xì)積分法; 衛(wèi)星編隊(duì)飛行; Runge-Kutta方法

北京大學(xué)學(xué)報(bào)（自然科學(xué)版）第52卷第4期2016年7月

Acta Scientiarum Naturalium Universitatis Pekinensis, Vol. 52, No. 4 (July 2016)

衛(wèi)星編隊(duì)飛行指若干個(gè)小衛(wèi)星按特定的編隊(duì)飛行, 各衛(wèi)星協(xié)同工作, 形成一個(gè)分布式的編隊(duì)衛(wèi)星群。與大衛(wèi)星相比, 編隊(duì)衛(wèi)星群成本低、可靠性強(qiáng)、靈活性高、編隊(duì)自由, 能完成大衛(wèi)星所不能完成的特殊任務(wù)。因此, 編隊(duì)衛(wèi)星技術(shù)在導(dǎo)航、通信、觀測(cè)、偵查等軍事和民用方面都扮演著重要角色, 具有廣闊的應(yīng)用前景[1]。

編隊(duì)的控制是編隊(duì)衛(wèi)星群運(yùn)行的重要前提, 包括編隊(duì)調(diào)整與編隊(duì)保持, 編隊(duì)控制的基礎(chǔ)則是衛(wèi)星間的相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)關(guān)系。1960 年, Clohessy 等[2]最早開(kāi)展衛(wèi)星飛行的相對(duì)運(yùn)動(dòng)動(dòng)力學(xué)關(guān)系研究, 并針對(duì)空間兩臨近飛行器交會(huì)問(wèn)題, 提出 Clohessy-Wiltshire方程(C-W方程)。C-W方程是線(xiàn)性化、常系數(shù)的微分方程組, 形式簡(jiǎn)單, 存在解析解[1], 在早期的衛(wèi)星編隊(duì)飛行和衛(wèi)星空間交會(huì)的動(dòng)力學(xué)與控制中得到廣泛應(yīng)用[3-4]。然而, 萬(wàn)有引力與衛(wèi)星軌道半徑之間的關(guān)系是非線(xiàn)性的, 因此描述衛(wèi)星編隊(duì)飛行的動(dòng)力學(xué)方程是非線(xiàn)性方程[5]。相對(duì)于衛(wèi)星的軌道半徑, 編隊(duì)衛(wèi)星群的衛(wèi)星間距很小, 若將萬(wàn)有引力引起的非線(xiàn)性項(xiàng)做泰勒展開(kāi), 高階項(xiàng)部分?jǐn)?shù)量級(jí)相對(duì)較小, 其線(xiàn)性部分占有主導(dǎo)地位, 因此, 衛(wèi)星編隊(duì)飛行的動(dòng)力學(xué)方程呈現(xiàn)弱非線(xiàn)性的特點(diǎn)。在研究衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)時(shí), 若采用 C-W 模型而忽略非線(xiàn)性項(xiàng)的影響, 勢(shì)必帶來(lái)較大的誤差, 其后果是消耗更多的燃料以維持衛(wèi)星編隊(duì), 影響功能的充分發(fā)揮。對(duì)衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)方程的積分, 應(yīng)當(dāng)充分利用其弱非線(xiàn)性的特點(diǎn), 首先將其做泰勒展開(kāi)后寫(xiě)成半線(xiàn)性的形式(即方程包含線(xiàn)性部分和非線(xiàn)性部分), 然后用精細(xì)積分法[6-7]求解線(xiàn)性部分,用合適的積分方法求解高階項(xiàng)的非線(xiàn)性部分。目的是精確計(jì)算數(shù)量級(jí)較大的線(xiàn)性部分, 而由于高階非線(xiàn)性部分的數(shù)量級(jí)較小, 誤差也比較小。

針對(duì)線(xiàn)性常系數(shù)微分方程, 經(jīng)典常微分方程理論的主要困難是指數(shù)矩陣的計(jì)算。1994年, 鐘萬(wàn)勰等[6-7]針對(duì)指數(shù)矩陣的計(jì)算提出精細(xì)積分法(precise integration method, PIM), 用矩陣的加法代替矩陣乘法, 巧妙地避免了舍入誤差。精細(xì)積分法有多方面的優(yōu)勢(shì), 如可以采用大步長(zhǎng), 計(jì)算精度高, 穩(wěn)定性好[8]。另外, Ashi 等[9]比較了 6 種常用的指數(shù)矩陣的求解方法, 其中矩陣分解(matrix decomposition)方法無(wú)論在大步長(zhǎng)還是小步長(zhǎng)都能精確計(jì)算指數(shù)矩陣, 但當(dāng)矩陣獨(dú)立特征向量個(gè)數(shù)小于維數(shù)時(shí), 此方法的應(yīng)用比較困難。精細(xì)積分法則不受此限制, 無(wú)論步長(zhǎng)大小都能精確計(jì)算指數(shù)矩陣, 是一種有效的計(jì)算方法。

針對(duì)半線(xiàn)性微分方程, Hersch[10]于 1958 年首次提出指數(shù)積分的思想。Certaine[11]于 1960 年用常數(shù)變異公式, 首次給出指數(shù) Adams-Moulton 方法。Hochbruck 等[12]于 2010 年綜述了指數(shù)積分方法的發(fā)展, 并系統(tǒng)地總結(jié)了指數(shù)積分法的計(jì)算格式。由于指數(shù)積分方法的思想是精確地求解微分方程的線(xiàn)性部分, 而方程的剛性和高振蕩性通常表現(xiàn)在其線(xiàn)性部分, 所以在指數(shù)積分方法提出的初期, 主要用于解決剛性問(wèn)題或高振蕩問(wèn)題[13]。在指數(shù)積分方法的計(jì)算格式中, 指數(shù)矩陣往往難以精確計(jì)算, 尤其在維數(shù)較大的情況下。這限制了指數(shù)積分方法的進(jìn)一步發(fā)展, 因此在指數(shù)積分方法提出后的很長(zhǎng)時(shí)間都沒(méi)有得到重視。21世紀(jì)初, 由于新的針對(duì)指數(shù)矩陣的計(jì)算方法的提出[9], 指數(shù)積分方法得到快速發(fā)展。目前指數(shù)積分方法在各個(gè)領(lǐng)域得到廣泛應(yīng)用, 如延時(shí)問(wèn)題[14]、高振蕩問(wèn)題[15]、拋物型方程的求解[16]、李群算法[17]以及 Hamilton 系統(tǒng)的保辛算法的構(gòu)造[18]等。針對(duì)此方法的研究重點(diǎn)由計(jì)算格式的構(gòu)造逐漸轉(zhuǎn)移到其計(jì)算性能的研究, 如穩(wěn)定性[19-20]、收斂性[21]等。

本文在精細(xì)積分法與指數(shù)積分方法的基礎(chǔ)上,提出精細(xì)指數(shù)積分方法, 并將其用于弱非線(xiàn)性方程的求解, 為弱非線(xiàn)性方程的求解提供一種有效的計(jì)算方法, 為衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)仿真提供一種精確的算法。

1　衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)方程

衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)主要研究伴隨衛(wèi)星(follower satellite)相對(duì)于參考衛(wèi)星(reference satellite)的相對(duì)運(yùn)動(dòng)與受力之間的關(guān)系。如圖 1 所示, 假設(shè)參考衛(wèi)星在圓形軌道上運(yùn)行, 軌道半徑為 r = |r|, O2X2Y2Z2為軌道坐標(biāo)系, 其中O2Y2軸的方向與地心到參考衛(wèi)星的方向一致, O2X2方向?yàn)閰⒖夹l(wèi)星運(yùn)動(dòng)的反方向, O2Z2方向由右手螺旋定則確定。在軌道坐標(biāo)系中, 伴隨衛(wèi)星的位置矢量為r2= [x, y, z]T, 在不考慮攝動(dòng)力和控制力的情況下, 伴隨衛(wèi)星相對(duì)于參考衛(wèi)星運(yùn)動(dòng)的Lagrange函數(shù)為

其中,μ為地球引力常數(shù), m為伴隨衛(wèi)星的質(zhì)量, θ˙為參考衛(wèi)星的定常軌道角速度, 等號(hào)右邊第一項(xiàng)為伴隨衛(wèi)星的動(dòng)能, 第二項(xiàng)為勢(shì)能。由Lagrange方程可導(dǎo)出伴隨衛(wèi)星相對(duì)于參考衛(wèi)星的動(dòng)力學(xué)方程:

方程(2)的非線(xiàn)性項(xiàng)在于萬(wàn)有引力引起的方程右端的分母。衛(wèi)星編隊(duì)飛行的距離通常比參考衛(wèi)星的軌道半徑小得多, 即, 故方程(2)右端的分母部分表現(xiàn)為弱非線(xiàn)性。為了簡(jiǎn)化方程, 尋求解析解, 一般可將其做泰勒展開(kāi)并保留線(xiàn)性項(xiàng), 可得

由于方程的線(xiàn)性化會(huì)引入誤差, 當(dāng)誤差積累到一定程度時(shí)需要對(duì)衛(wèi)星施加控制, 從而消耗更多的燃料以維持編隊(duì), 甚至影響編隊(duì)衛(wèi)星群功能的實(shí)現(xiàn)。在求解衛(wèi)星編隊(duì)飛行的非線(xiàn)性方程時(shí), 為了方便精細(xì)指數(shù)積分方法的使用, 將方程(2)寫(xiě)成半線(xiàn)性的形式:

方程(4)右端的第一項(xiàng)為方程(2)的線(xiàn)性部分; 第二項(xiàng)和第三項(xiàng)合起來(lái)為方程(2)的高階非線(xiàn)性部分, 由于方程表現(xiàn)為弱非線(xiàn)性, 該部分的數(shù)值比線(xiàn)性部分小若干個(gè)數(shù)量級(jí)。

2　精細(xì)積分法

對(duì)于衛(wèi)星編隊(duì)飛行的線(xiàn)性化動(dòng)力學(xué)方程(式(3)), 可以利用精細(xì)積分法, 計(jì)算得到伴隨衛(wèi)星相對(duì)于參考衛(wèi)星的近似運(yùn)動(dòng)規(guī)律。精細(xì)積分法的詳細(xì)推導(dǎo)過(guò)程[6-7]如下。

將方程(3)降階成一階微分方程組:

其中,

若能精確求解矩陣 T, 則可以精確求解常系數(shù)齊次線(xiàn)性微分方程(5)。精細(xì)積分法的核心正在于矩陣T的精確求解, 求解思路是利用指數(shù)函數(shù)的加法原理, 將時(shí)間步長(zhǎng)τ分成2n等分, 令則有

由于tΔ是一個(gè)很小的時(shí)間段, 所以在 exp(LΔt)的泰勒展開(kāi)式中直接取其前 5 項(xiàng), 即可達(dá)到足夠的精度, 即

其中, I為單位矩陣。與I相比, T1是一個(gè)很小的矩陣。為防止計(jì)算過(guò)程中“大數(shù)吃小數(shù)”, 要避免T1與I直接相加。由方程(7)和(8)可知:

若有矩陣序列:

則根據(jù)式(9)有如下遞推公式:

通過(guò)求此矩陣序列得到 Tn+1, 最后與單位矩陣I相加, 便可計(jì)算矩陣 T 的值。在此過(guò)程中, 巧妙地運(yùn)用指數(shù)函數(shù)的性質(zhì), 避免了嚴(yán)重的舍入誤差, 所以這種算法有很高的精度, 一般取n=20可保證算法良好的性能。由于計(jì)算過(guò)程中在每個(gè)時(shí)間步長(zhǎng)τ 內(nèi)插入220個(gè)點(diǎn), 所以即使步長(zhǎng)比較大, 也不影響算法的精確性。精細(xì)積分法求解一般的常系數(shù)齊次線(xiàn)性微分方程, 能夠達(dá)到計(jì)算機(jī)意義上的精度。

3　精細(xì)指數(shù)積分方法

精細(xì)積分法只能處理衛(wèi)星編隊(duì)飛行的線(xiàn)性化動(dòng)力學(xué)方程。為了更加精確地求解伴隨衛(wèi)星相對(duì)于參考衛(wèi)星的動(dòng)力學(xué)規(guī)律, 還需考慮其非線(xiàn)性因素。將半線(xiàn)性微分方程(4)降階, 且寫(xiě)成如下形式:

其中, L與方程(5)中的L相同, N(u, t)為方程的非線(xiàn)性部分, 當(dāng)方程為弱非線(xiàn)性時(shí)有

對(duì)時(shí)間離散后, 方程(11)的精確解可寫(xiě)為如下形式:

其中, 右端的積分項(xiàng)稱(chēng)為Duhamel積分。用不同的方法求解, 可得到不同的指數(shù)積分方法, 目前應(yīng)用較廣的是指數(shù) Runge-Kutta 方法。s 級(jí)的指數(shù) Runge-Kutta方法格式如下:

其中, aij(Lτ)和bi(Lτ)均為待定矩陣, 且為 L 和τ 的函數(shù), ci為待定常數(shù)。當(dāng)N(u, t)=0時(shí), 指數(shù)Runge-Kutta 方法退化為常系數(shù)齊次線(xiàn)性微分方程的精確解, 當(dāng) L=0 時(shí), 此方法退化為經(jīng)典的 Runge-Kutta方法。

根據(jù) aij(Lτ)和 bi(Lτ)給出的方法不同, 指數(shù)Runge-Kutta方法可以分成不同的類(lèi)型, 其中包括兩種最常用的類(lèi)型: 一種是積分因子法(integration factor, IF), 也稱(chēng)為 Lawson 方法[22], 這種方法的待定矩陣完全根據(jù)已有的 Runge-Kutta 方法確定, 構(gòu)造和應(yīng)用都比較簡(jiǎn)單; 另一種是指數(shù)時(shí)間微分方法(exponential time differencing, ETD), 此方法通過(guò)常數(shù)變異公式, 將非線(xiàn)性部分用代數(shù)插值多項(xiàng)式代替, 需要定義φ函數(shù)[12,19]。本文只對(duì)積分因子法做簡(jiǎn)單介紹。

積分因子法中的 aij(Lτ)和 bi(Lτ)按如下規(guī)律[19]給出:

員工是企業(yè)發(fā)展的關(guān)鍵，企業(yè)發(fā)展的動(dòng)力是廣大員工，企業(yè)員工的管理工作十分重要。當(dāng)然，在知識(shí)經(jīng)濟(jì)的大環(huán)境下，對(duì)員工管理所存在的各種問(wèn)題，人們逐漸有了全新的認(rèn)識(shí)。當(dāng)前社會(huì)，很多企業(yè)都十分關(guān)心企業(yè)客戶(hù)的滿(mǎn)意度。隨著整個(gè)市場(chǎng)競(jìng)爭(zhēng)活動(dòng)的不斷加劇，越來(lái)越多的企業(yè)認(rèn)識(shí)到完善的經(jīng)營(yíng)理念并非只有外部因素。因此，企業(yè)經(jīng)營(yíng)者不僅需要關(guān)注外部環(huán)境的要求，還需要關(guān)注客戶(hù)與企業(yè)內(nèi)部的員工管理活動(dòng)，探索適合當(dāng)前企業(yè)經(jīng)營(yíng)的新道路。

其中, aij, bi, ci與經(jīng)典Runge-Kutta方法的系數(shù)相同,即給定一種經(jīng)典的Runge-Kutta方法就有一種積分因子法與之對(duì)應(yīng), 且當(dāng) N(u, t)=0 時(shí)積分因子法退化為相應(yīng)的經(jīng)典的Runge-Kutta方法。例如, 經(jīng)典四級(jí)四階的Runge-Kutta方法所對(duì)應(yīng)的積分因子法如下:

4　衛(wèi)星編隊(duì)飛行數(shù)值仿真

在本節(jié)中, 先通過(guò)一個(gè)弱非線(xiàn)性的算例, 驗(yàn)證本文提出的積分因子法的有效性, 然后將此方法用于衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)方程的求解。

4.1積分因子法有效性的驗(yàn)證

考慮以下弱非線(xiàn)性微分方程:

方程有解析解:

當(dāng) a 取值較小時(shí), 方程的非線(xiàn)性部分?jǐn)?shù)量級(jí)比較小, 此處取 a=0.1。分別采用經(jīng)典的四級(jí)四階Runge-Kutta 方法(RK4_4)、八級(jí)六階 Runge-Kutta方法(RK8_6)[23]、三級(jí)六階隱式 Runge-Kutta 方法(RK3_6)[23]、與 RK4_4 對(duì)應(yīng)的積分因子法(IF4_4)、與 RK3_6 對(duì)應(yīng)的積分因子法(IF3_6), 在不同的時(shí)間步長(zhǎng)下求解方程(16), 結(jié)果如圖2所示。

圖2為以上各種方法相對(duì)誤差絕對(duì)值的最大值與積分時(shí)間步長(zhǎng)之間的關(guān)系。由圖 2 可知, RK4_4與IF4_4、RK3_6與IF3_6基本上平行, 說(shuō)明積分因子法與對(duì)應(yīng)的 Runge-Kutta 方法具有相同的階數(shù)。在相同的步長(zhǎng)時(shí), IF4_4比RK4_4、IF3_6比RK3_6的計(jì)算誤差都低4個(gè)數(shù)量級(jí)左右, 說(shuō)明在弱非線(xiàn)性情況下, 指數(shù) Runge-Kutta 方法比 Runge-Kutta 方法在精度上有比較明顯的優(yōu)勢(shì)。將 IF4_4與 RK3_6 進(jìn)行比較, 在較大步長(zhǎng)時(shí), IF4_4 方法比RK3_6 更精確; 在中等步長(zhǎng)時(shí), 兩種方法精度相當(dāng);而在較小步長(zhǎng)時(shí), RK3_6更有優(yōu)勢(shì)。

為了分析方程非線(xiàn)性部分的數(shù)值大小與各種方法求解相對(duì)誤差之間的關(guān)系, 取步長(zhǎng)為0.1τ=, a在0.01~1 之間變化, 同樣用上述5種方法求解方程(16), 結(jié)果如圖3所示。

由圖 3 可知, 對(duì)于普通的 Runge-Kutta 方法,方程非線(xiàn)性的強(qiáng)弱對(duì)求解結(jié)果的相對(duì)誤差影響不大, 然而, 精細(xì)指數(shù)積分方法求解結(jié)果的相對(duì)誤差隨非線(xiàn)性部分的數(shù)值增大而增加。方程的非線(xiàn)性越弱, 精細(xì)指數(shù)積分方法求解結(jié)果越準(zhǔn)確。與同階的Runge-Kutta 方法相比, 精細(xì)指數(shù)積分方法在精度上有優(yōu)勢(shì)。

4.2衛(wèi)星編隊(duì)調(diào)整動(dòng)力學(xué)仿真

在衛(wèi)星編隊(duì)飛行的非線(xiàn)性動(dòng)力學(xué)方程中, 對(duì)于Runge-Kutta 方法, 方程(2)與方程(4)沒(méi)有本質(zhì)上的差別, 為了便于計(jì)算, 可直接積分方程(2); 對(duì)于指數(shù)Runge-Kutta方法, 先將方程寫(xiě)成半線(xiàn)性形式(即方程(4)), 再進(jìn)行計(jì)算。由于非線(xiàn)性動(dòng)力學(xué)方程沒(méi)有解析解, 所以用RK3_6和RK8_6兩種高階算法的計(jì)算結(jié)果作為參考, 主要比較RK4_4與IF4_4的計(jì)算結(jié)果。另外, 用精細(xì)積分法(PIM)求解線(xiàn)性化的動(dòng)力學(xué)方程(方程(3)), 研究方程的線(xiàn)性化帶來(lái)的誤差。

假設(shè)在一次編隊(duì)調(diào)整過(guò)程中, 參考衛(wèi)星的軌道半徑為8000 km, 伴隨衛(wèi)星相對(duì)于參考衛(wèi)星的初始狀態(tài)為

積分步長(zhǎng)為τ =500 s, 經(jīng)過(guò) 4 小時(shí)的運(yùn)行, 軌跡圖如圖4所示。的結(jié)果, 與其他曲線(xiàn)均有明顯差距。由圖 4 可知,經(jīng)過(guò) 4 小時(shí)運(yùn)行后, 線(xiàn)性化方程的計(jì)算誤差為 102m 量級(jí)。從圖 4 的局部放大圖中可看到, IF4_4 的計(jì)算軌跡與RK3_6和RK8_6的計(jì)算軌跡基本上重合, 而RK4_4的誤差為 10 m 量級(jí)。對(duì)比IF4_4與RK4_4 的計(jì)算結(jié)果可知, 在衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)方程的求解中, 精細(xì)指數(shù)積分方法能提供更好的計(jì)算精度。

圖 4 中PIM即為精細(xì)積分法求解線(xiàn)性化方程

4.3衛(wèi)星編隊(duì)保持動(dòng)力學(xué)仿真

當(dāng)參考衛(wèi)星的軌道半徑為 8000 km 時(shí), 伴隨衛(wèi)星相對(duì)于參考衛(wèi)星的初始狀態(tài)為即可使伴隨衛(wèi)星相對(duì)于參考衛(wèi)星的運(yùn)動(dòng)軌道為周期軌道。分別用RK3_6, RK4_4和IF4_4方法求解動(dòng)力學(xué)方程, 取積分步長(zhǎng)為τ = 600 s, 積分時(shí)間為100小時(shí), 得到伴隨衛(wèi)星相對(duì)于參考衛(wèi)星的運(yùn)動(dòng)軌跡圖(圖5)。

由圖 5 可知, 六階的 RK3_6 方法和四階的IF4_4方法能精確地模擬周期軌道, RK4_4方法求解的周期軌道則出現(xiàn)較大的誤差。由此可知, 在弱非線(xiàn)性問(wèn)題的求解上, 與同階的 Runge-Kutta 方法對(duì)比, 指數(shù)積分方法能提供較高的精度, 因此在求解時(shí)可以采用較大的步長(zhǎng)。

5　結(jié)論

本文將精細(xì)積分法用于指數(shù)積分方法中的指數(shù)矩陣的求解, 得到精細(xì)指數(shù)積分方法, 并將其用于弱非線(xiàn)性方程的求解。針對(duì)弱非線(xiàn)性問(wèn)題, 例如衛(wèi)星編隊(duì)飛行的動(dòng)力學(xué)方程, 精細(xì)積分法不再適用。如果忽略其弱非線(xiàn)性的性質(zhì), 直接用 Runge-Kutta方法求解, 則線(xiàn)性部分也會(huì)引入誤差。

精細(xì)指數(shù)積分方法能精確求解方程的線(xiàn)性部分, 并用合適的算法求解非線(xiàn)性部分, 保證了算法的精度。與同階的 Runge-Kutta 方法對(duì)比, 方程的非線(xiàn)性越弱, 精細(xì)指數(shù)積分方法在精度上越有優(yōu)勢(shì)。在衛(wèi)星編隊(duì)飛行的動(dòng)力學(xué)仿真中, 精細(xì)指數(shù)積分方法取得了比同階的Runge-Kutta方法更精確的結(jié)果。

[1] 李俊峰, 高云峰, 寶音賀西, 等. 衛(wèi)星編隊(duì)飛行動(dòng)力學(xué)與控制研究. 力學(xué)與實(shí)踐, 2002, 24(2): 1-6

[2] Clohessy W H, Wiltshire R S. Terminal guidance system for satellite Rendezvous. Journal of the Aerospace Sciences, 1960, 27(9): 653-658

[3] 周建平. 天宮一號(hào)/神舟八號(hào)交會(huì)對(duì)接任務(wù)總體評(píng)述. 載人航天, 2012, 18(1): 1-5

[4] 王博, 鄧子辰, 李文成, 等. Hill’s方程的Magnus積分. 西北工業(yè)大學(xué)學(xué)報(bào), 2011, 29(6): 988-991

[5] 林來(lái)興. 交會(huì)對(duì)接動(dòng)力學(xué)模型和動(dòng)力學(xué)特性. 中國(guó)空間科學(xué)技術(shù), 1994, 14(3): 47-53

[6] 鐘萬(wàn)勰. 結(jié)構(gòu)動(dòng)力方程的精細(xì)時(shí)程積分法. 大連理工大學(xué)學(xué)報(bào), 1994, 34(2): 131-136

[7] Zhong Wanxie, Williams F W. Precise time step integration method. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, 1994, 208(6): 427-430

[8] Zhong Wanxie. On precise integration method. Journal of Computational and Applied Mathematics, 2004, 163(1): 59-78

[9] Ashi H A, Cummings L J, Matthews P C. Comparison of methods for evaluating functions of a matrix exponential. Applied Numerical Mathematics,2009, 59(3): 468-486

[10] Hersch J. Contribution à la méthode des équations aux différences. Zeitschrift für angewandte Mathematik und Physik ZAMP, 1958, 9(2): 129-180

[11] Certaine J. The solution of ordinary differential equations with large time constants. Mathematical Methods for Digital Computers, 1960, 1: 128-132

[12] Hochbruck M, Ostermann A. Exponential integrators. Acta Numerica, 2010, 19(1): 209-286

[13] Rahrovani S, Abrahamsson T, Modin K. An efficient exponential integrator for large nonlinear stiff systems Part 1: Theoretical investigation.Conference Proceedings of the Society for Experimental Mechanics Series, 2014, 2: 259-268

[14] Xu Yang, Zhao Jingjun, Sui Zhenan. Exponential Runge-Kutta methods for delay differential equations. Mathematics and Computers in Simulation, 2010, 80(12): 2350-2361

[15] Frénod E, Hirstoaga S A, Sonnendrücker E. An exponential integrator for a highly oscillatory Vlasov equation. Discrete and Continuous Dynamical Systems: Ser S, 2015, 8(1): 169-183

[16] Luan V T, Ostermann A. Explicit exponential Runge-Kutta methods of high order for parabolic problems. Journal of Computational and Applied Mathematics, 2014, 256: 168-179

[17] Celledoni E, Marthinsen A, Owren B. Commutatorfree Lie group methods. Future Generation Computer Systems, 2003, 19(3): 341-352

[18] Miyatake Y. An energy-preserving exponentiallyfitted continuous stage Runge–Kutta method for Hamiltonian systems. BIT Numerical Mathematics, 2014, 54(3): 777-799

[19] Maset S, Zennaro M. Unconditional stability of explicit exponential Runge-Kutta methods for semilinear ordinary differential equations. Mathematics of Computation, 2009, 78: 957-967

[20] Maset S, Zennaro M. Stability properties of explicit exponential Runge-Kutta methods. IMA Journal of Numerical Analysis, 2013, 33(1): 111-135

[21] Berland H, Owren B, Skaflestad B. B-Series and order conditions for exponential integrators. SIAM Journal on Numerical Analysis, 2005, 43(4): 1715-1727

[22] Lawson J D. Generalized Runge-Kutta processes for stable systems with large Lipschitz constants. SIAM J Numer Anal, 1967, 4: 372-380

[23] 馬振華. 現(xiàn)代應(yīng)用數(shù)學(xué)手冊(cè): 計(jì)算與數(shù)值分析卷.北京: 清華大學(xué)出版社, 2005

Precise Exponential Integrator and Its Application in Dynamics of Spacecraft Formation Flying

DENG Zichen?， LI Qingjun

Department of Engineering Mechanics, Northwestern Polytechnical University, Xi’an 710072; ? E-mail: dweifan@nwpu.edu.cn

The dynamic equations of spacecraft formation flying are weakly nonlinear equations since the distance between spacecrafts is quite small compared with the orbital radius of the spacecrafts. To solve weakly nonlinear equations effectively, a precise exponential integrator (PEI) was proposed. Precise integration method (PIM) was applied to calculate exponential function in the formulas of exponential integrators (EI). Firstly, PEI was validated by solving a weakly nonlinear equation compared with Runge-Kutta method. Secondly, the dynamic equations of spacecraft formation flying were obtained through Lagrange equations, and then the equations were tansfered into semi-linear form. Ultimately, PEI and Runge-Kutta method were comparatively used to solve these equations. Through numerical analysis, PEI gave higher precision of the dynamic equations of spacecraft formation flying, indicating that PEI can be applied to other weakly nonlinear problems as well.

exponential integrator; precise integration method; spacecraft formation flying; Runge-Kutta method

O316

10.13209/j.0479-8023.2016.069

國(guó)家自然科學(xué)基金(11432010)資助

2015-10-07;

2016-02-03; 網(wǎng)絡(luò)出版日期: 2016-07-14