Cowell數值積分器的變步長與自起步方法

夏 炎, 田 波

(1.銅仁學院大數據學院，銅仁 554300；2.中國科學院紫金山天文臺，南京 210008)

數值積分是一個經典問題，自20世紀70年代起，伴隨著計算機技術的發展，對積分器的研究一直沒有停止過。近幾年仍不斷有新的積分方法提出，計算高效、程序實現簡單、適用性廣是總的發展趨勢[1-5]。積分器分為單步法和多步法，前者每步積分都需要計算多次被積函數，且步長較小；多步法積分器對被積函數的計算次數少，積分步長大，在計算效率上具有優勢。

多步法積分器又分為等間距和可變間距。可變間距積分器可以靈活控制積分步長，但每積分一步都需要重新計算積分器系數，程序實現較為復雜[6]。而等間距的多步法積分器程序實現相對簡單，使得它在人造衛星、自然天體的軌道計算中應用廣泛[7-11]。然而在步長變換上一直都是難點，限制了其在大偏心率軌道、天體密近交匯等情況下的使用。如果可以為這類積分器解決變步長的困難，將對實際應用有很大幫助。

多步法積分器還有一個缺點就是積分器的起步比較復雜，需要在初始點附近求解多個點的狀態后，才能啟動正常積分。在實際應用中廣泛使用的方法是求解各個步點的解析近似解然后迭代，或直接借助單步法積分器求解各步點再起步[6, 12-14]。如果可以擺脫對近似解或單步法積分器的依賴，將使積分器的使用更加方便和流暢。

自20世紀70年代起，紫金山天文臺行星室就一直使用等間距的Cowell多步法積分器進行相關領域的研究[15-17]，在此過程中也對積分器本身提出了一些改進措施，文獻[14]做過一次總結性的敘述。筆者在上述文獻的基礎上對Cowell積分器做了進一步改進。

首先推導了高階的非整數步點坐標速度的求解系數，使得高精度計算非整數步點狀態成為可能，以此為基礎實現了積分步長的變換操作。隨之而來的問題是變換步長的時機與步長大小的選取，提供了積分器截斷誤差的計算方法，可以以此為依據決定是否改變積分步長，以及預估步長改變之后的截斷誤差大小。利用二體問題模型對所給變步長方法做了檢驗，結果表明即便是幾十萬次的反復變步長操作也不會引入誤差。

其次實現了積分器的逐階迭代起步。逐階迭代起步是理論上成熟的方法，但在實際應用中很少使用。其主要原因是其中要用到積分器每一階的系數，而低階積分器使用較少，其積分器系數也不易查到。給出了各階積分器系數、算法流程與程序實例，方便參考使用。

最后介紹一種減小積分器舍入誤差的措施，可以大幅減小積分器誤差。有助于解決高精度的深空探測數據處理、超長期的行星系統演化問題。積分器精度優于以能量保持性占優的IAS15積分器[5]。

涉及的積分器系數以及相關算法的示例程序源代碼已經上傳到github網站[18-19]。

1 基本算法

需要求解的基本方程為

(1)

為了推導積分公式有以下算子：

移位算子Ef(t)=f(t+h)；

如文獻[4]中的推導，算子之間存在如下關系：

E=ehD。

于是有：

(2)

(3)

再根據中差算子的運算[14]就可以把坐標和速度的計算與已知步點的函數值聯系起來。

(4)

如果函數表和和分表已知，那就可以根據式(2)和式(3)計算各步點的坐標速度，最終坐標的計算簡化如下：

(5)

式(5)中：yk是第k個步點的坐標，k∈[-7,7]；U是一個8×13的常系數矩陣。速度的計算也有類似的公式，將在后文統一給出。

要進行一次積分運算，遵循如下流程：首先構建13個已知步點的函數值，使用較多的方法是用解析近似解，迭代至收斂；然后利用式(5)和初始坐標速度，以及式(4)，完成和分表；最后可以逐步外推，外推一步的流程如下。

(1)利用式(4)計算+7步點和分值。

(2)利用式(5)計算+7步點的坐標速度，進而計算函數值，稱為預報。

(3)函數表、和分表移動一個步點。

(4)利用式(5)計算+6步點的坐標速度，進而計算函數值，稱為修正。

2 變步長的實現

積分器改變步長就是重新構建13個已知點，并使它們之間的距離縮小或放大。首先給出非整數步點狀態的高精度求解，用來計算任意時刻的被積函數，以此構建13個新步點，從而實現步長變換的操作；然后給出截斷誤差的計算，作為步長選取和變換時機的判據。

2.1 非整數步點坐標速度的計算

以往對于非整數步點狀態的求解只是用來作為輸出工具，其結果不參與積分迭代，要求的精度不高。文獻[14]中推導了8階的積分器，但非整數步點的計算精度只準確到了4階。筆者針對12階的積分器推導同樣準確到12階非整數步點的系數。

設實數n∈(-1,1)，則在整數步點k附近的k+n處有：

fk+n=Enfk=enhDfk=

(6)

使用(hD)-2作用在fk上即得整數步點的坐標，現在把它作用在fk+n上就可以得到非整數步點坐標的計算公式,即

(hD)-2fk+n=[(hD)-2+n(hD)-1+

(7)

然后類似于整數步點的處理，可以得到使用函數表計算坐標的方法，其中需要15個8×13的常系數矩陣，記為U15,8,13；最終對于步點k∈[-7,7]附近的位置k+n處，非整數步點的坐標計算公式如下：

當k≥0時，

(8)

當k≤0時，

(9)

式中：

(10)

同樣的方法可以得到速度計算公式。

當k≥0時，

(11)

當k≤0時，

(12)

若n=0，則不需要再對j求和，只需要計算j取值最小的單項就可以了，也就回到了整數步點的計算公式。

2.2 變步長的操作

變步長的操作核心是按照新步長重新構建13個步點。在解決了非整數步點的坐標速度的求解問題之后，步長的縮小就是簡單的。步長減半的操作如下。

(1)保存中心步點的坐標速度。

(2)計算-2.5、-1.5、-0.5、0.5、1.5、2.5步點處的坐標速度，再計算右函數，結合已知的-3到+3的整數步點的加速度，就構成了13個已知點。

(3)使用保存的中心步點坐標速度，根據式(5)確定和分表。

為了便于擴大步長，對傳統的積分器使用方法做了一點改變，增加了保存的函數值的個數。傳統方法只要保存13個步點的函數值，就可以使積分器外推。現在改為保存25個步點，這樣就可以很方便地實現步長加倍的操作。而在今天的計算條件下，多保留的這些步點并不會造成計算負擔。

對上述變步長方法進行了測試，使用日地二體模型，對地球軌道積分1 000圈，積分結果與二體模型解析結果比較，如圖1所示。

圖1 變步長對積分精度的影響Fig.1 The influence of change step size on precision

從圖1中可以看出，變步長的操作不會引入額外誤差，而造成精度損失。采用固定積分步長1、2、4 d做積分運算。對變步長做了兩個測試：在1 d與2 d之間反復切換，和在2 d與4 d之間反復切換。步長為1 d和2 d之間切換時，增大步長365 249次，減半365 250次；步長為2 d和4 d之間切換時，增大步長182 624次，減半182 625次。從圖1中可以看出，雖然有幾十萬次的變步長操作，但并未對計算精度造成影響。(計算中使用的變量類型是IEEE規范的雙精度浮點數類型。)

2.3 變步長的判據

解決了積分器變步長的具體操作后，另一個問題是變步長的時機判斷，在積分過程中怎樣判斷當前的積分步長是否合適，需要增大還是減小。

文獻[6]中討論了Gauss-Jackson積分器改變步長的判據，是使用預報值與改正值的差作為判據，理由是預報值比修正值少使用了一個步點，所以可以認為預報值精度低1階。但Gauss-Jackson積分器是基于后差算子的，而在由中差算子導出的Cowell積分器中，沒有明顯的證據說預報值比修正值少使用了一個步點。所以使用預報值與修正值的差作為判據就有些牽強。

利用式(2)可以很嚴格地給出計算積分的12階截斷誤差為5.104×10-7δ12，結合文獻[4]對差分的計算，最終可以得到12階截斷誤差的計算方法如下：

(13)

式(13)中：A是常系數數組，其值為{1，-12，66，-220，495，-792，924，-792，495，-220，66，-12，1}。

注意到式(13)中截斷誤差的計算只和函數表有關，而在上一節中，把函數表保留的步點增加了一倍，使得擴大步長的操作變得簡便。在這里，這些多保留的步點也就可以用來預判擴大步長后的截斷誤差，即可以事先計算積分步長增大一倍后的截斷誤差會是多大，這有助于實際使用。

3 積分器的自起步

積分器要解決的是給定初始狀態的常微分方程，已知狀態只有一個時間點。但是要讓多步法積分器開始工作，需要首先構建一系列的已知點，計算這些點的函數值，如12階的積分器就需要13個這樣的已知點，而為了求解這13個點的函數值，通常是借助積分器以外的工具。

第一種方法是求助于單步法積分器，該方法應用比較普遍[20]。另一種方法是為各個步點尋求近似解析解，然后迭代至給定精度。關于解析近似解的來源，有的采用二體模型下的解[12, 14]，有的使用級數展開方式考慮一定的攝動因素[6]。

無論是借助單步法積分器，還是借助解析近似解都讓積分器欠缺單獨工作的能力，在文獻[21]中有針對2階積分器的起步方法研究，但并不能推廣到更高階的情況。在此給出一種逐階迭代的方式來實現積分器的自起步，具體流程如下。

(3)把得到的3個點的函數值，作為2階積分器的初始近似，迭代至收斂；然后向前向后各外推一個點，作為4階積分器的初始近似，再迭代至收斂。

(4)類似前一步，從4階到6階，再到8階……一直到12階。

算法的實現需要用到由2階至12階，每隔兩階的所有積分器系數，系數推導方法如前文所述，具體系數矩陣可在示例源碼中找到。

4 舍入誤差的減小

如何減小積分器誤差是個永恒的問題[22]。Al-Mohy等[23]討論過浮點數求和的舍入誤差問題，并給出了幾種減小舍入誤差的方法；Grazier等[24]把其中方法應用在了積分器上。紫金山天文臺張家祥等長期從事數值計算研究[25]，針對Cowell積分器也提出過一種減小舍入誤差的措施，簡單有效，在科研工作中一直被使用，文獻[14]中曾對此做過說明。筆者為了使闡述的積分器更加完善，也做一下介紹。

浮點數的大數和小數相加減會丟失小數的精度，特別的，積分器的外推是依照式(4)進行的，而其中的和分值″F和函數值F通常都有幾個數量級的相差，它們的加減是積分器誤差的主要來源，這個誤差會引入迭代中去，而且每一步迭代又都有新的誤差產生。

使用兩個浮點數變量存儲一個和分值的方法，以增加和分值存儲的有效數字位數。具體的操作是讓和分值乘以一個系數(記為Z)，使乘積的整數位具有一定的有效位數，然后把該乘積的整數部分和小數部分分別用兩個變量存儲。當有數值要與和分值做加法運算時，把該數值也乘以Z，然后整數與整數部分相加，小數與小數部分相加，并把小數部分產生的整數清除，加到整數部分里。

上述方法可以大幅提高積分器精度。對外太陽系天體做長期積分，考察其能量保持性，積分108個木星周期，能量差在10-13的量級，這個結果比以能量保持性占優的IAS15積分器[2]還要好。但是能量保持的精度并不代表積分器的精度[26]，所以仍然使用日地二體模型給出一個對比結果，如圖2所示。

圖2中縱軸是積分結果和解析解的位置誤差，橫軸是積分時間。可以看出增加字長措施大大減小了積分誤差，使得積分109圈以后地球仍然保持了它大體的位置，精度在10-3。計算中使用的變量類型是IEEE規范的雙精度浮點數，積分步長為1 d。

圖2 增加和分存儲字長的效果Fig.2 The effect of extending storage length of summations

5 結論

Cowell積分器簡單易用、計算效率高，在天體力學中有著廣泛的應用，但也存在著變步長困難、積分起步相對復雜的缺點。針對這些問題做了改進，并介紹了一些使積分器更加完善的措施。主要改進有以下幾點。

(1)提供了高精度計算非整數步點狀態的具體方法。如果沒有非整數步點的計算，積分器在積分過程中就是一些離散的點，有了非整數步點的狀態的計算，使得積分器更像是在時間軸上行進的一段連續函數。這對于處理某些問題，諸如光行時修正等，是非常便利的。

(2)提供了一種變換步長的方法，這種變步長方法不會引入額外的積分誤差，并介紹了變步長的判據。積分器具備了變步長能力后，將大大拓展它的使用空間，例如深空探測中探測器飛掠大行星的問題，如果采用固定步長，將不得不整條軌道都采取很小的積分步長，徒增計算量。

(3)改進了積分器的起步流程，使得積分器可以單獨工作，不需要借助于解析近似解或其他積分器起步。

(4)介紹了一種增加和分值存儲字長的方法，該方法可以大幅減小舍入誤差。使得二體模型下的被積天體在計算109圈后，仍然保持可信賴的位置精度。

對Cowell積分器的改進可為其他積分器提供參考。