天體力學第10課星曆表1問題5：克普勒方程

pang

天體力學第10課星曆表1問題5：克普勒方程

這問題分為兩部份，(1)克普勒方程的推導 及 (2)克普勒方程的解法。

克普勒方程的推導在「牧夫論壇」曾討論過，他們提出兩方法，方法一由wenlianyi君提出，較為複雜，整理後我再介紹給同好；方法二由zhjj 君提出，和我在下列附件中提出的方法相同。

克普勒方程的解法在「牧夫論壇」提及的課本「天文算法」中有詳細討論，此書有英文本，亦有中文翻譯本，由许剑伟译著，同好可在網站 http://www.uushare.com/user/mglychd/files/1191390下載。

若同好發現克普勒方程有其他推導或解法，請與我們分享。
謝謝！
彭祿勝  謹啟   2009-11-22
-----------------------------------------------------------------------------------------
附件：克普勒方程(Kepler equation)  M=E-e˙sinE 的推導

附圖中，設
E：偏近點角∠QOW (Eccentric anomaly)
M：平近點角(Mean anomaly)
M0：曆元平近點角(Mean anomaly of the epoch)
t：時間
τ：過近日點時間或曆元(Epoch)
T：週期(Period)
n：平均角速度(Mean angular velocity)
又 M = M0+nt = n(t-τ) = 2π(t-τ)/T ; n=2π/T

由幾何特性得：
SWP面積
= (b/a)˙SWQ
= (b/a)˙(扇形OWQ-ΔOSQ)
= (b/a)˙(0.5˙a^2˙E - 0.5˙a^2˙e˙sinE)
= 0.5ab˙(E-e˙sinE)

由克普勒第二定律： SWP面積/(πab) = (t-τ)/T
代入得： 0.5ab˙(E-e˙sinE) /(πab) = (t-τ)/T
即 E-e˙sinE = M
-----------------------------------------------------------------------------------------

pang

克普勒方程的解法： E = M + e*sin(E)

展開法：(参見「天体力学引论」第84頁，易照华著，科学出版社，1978年第1版)。
E = M + e*sin(M) + 0.5*e^2*sin(2M) + ...

「天文算法」(Astronomical Algorithms, Meeus)中譯本第29章提出4種方法，(1)迭代法、(2)牛頓公式法、(3)二分法和 (4)近似法。

(1) 迭代法：公式如下
E0 = M
E1 = M + e*sin(E0)
E2 = M + e*sin(E1)
. . . . . . . . . .

(2) 牛頓公式法：公式如下
E0 = M
E1 = E0 - (E0-M-e*sin(E0))/(1-e*cos(E0))
E2 = E1 - (E1-M-e*sin(E1))/(1-e*cos(E1))
. . . . . . . . . . . . . . . . . . . . .

迭代法收斂慢，牛頓公式法收斂快；當e接近1時，迭代法收斂極慢，牛頓公式法亦有些因難。

本人有這樣意見，太陽系天體的橢圓軌道偏心率e值極少接近1 (e值極接近1時已被算為拋物線軌道)，可一律用牛頓公式解克普勒方程，但在每一循環中加入以下課語句，令E值在0至2*PI之間(PI為圓周率) 。
E = MOD(E,2*PI)
IF (E<0)  E = E + 2*PI
MOD為餘數函數，例 MOD(7,3)=1，即7除以3餘數為1。

參照「天文算法」第164-165頁中，若e=0.999，M=7˚度，牛頓法不調整下需47步才得E=52.2702615˚，我利用牛頓法餘數調整後袛需7步亦得相同結果；若e=0.99999，M=7˚，調整後牛頓法亦是用了7步得 E=52.385629080494

彭祿勝  謹啟   2009-11-29

LLH

可以藉設定不同的初始值以減少迭代次數。

例如可試
a = sqrt(8*(1-e)/e)
u = arc sinh(3*M/(a*(1-e)))
那麼選 E0 = a*sinh(u/3) 代替 E0 = M
無論e是否接近1，均可令迭代次數大減。  (M 須以弳制而非度分秒表示， sinh 及 arc sinh 分別為雙曲正弦及反雙曲正弦函數，在 Excel 試算表中也有提供)

mca

1985 年 S&T 雜誌 有一  BASIC Language Program 計算 Kepler Equation  E = M + e *sin (E)
http://media.skyandtelescope.com/binary/kepler.bas

我沒有細看 program (放下 BASIC 太久了，有點忘記所學)，或者這 BASIC 可減低人力計算的繁複步驟，1985 是 PC (Intel 386 CPU) 發育年代，現在的 PC 更勁百倍，可試試以 PC 解。

我又翻查「天文演算法」 原著 Astronomical Algorithms Chapter 30，它也介紹了四個計算方法，其中的二分法 (binary search method, p.206) 與 S&T 的 BASIC 相似。

Alan