天體力學第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
-----------------------------------------------------------------------------------------
天體力學第10課星曆表1問題5:克普勒方程
天體力學第10課星曆表1問題5:克普勒方程
- 附加檔案
-
- ZZ.GIF (12.24 KiB) 已瀏覽 5665 次
克普勒方程的解法: 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
展開法:(参見「天体力学引论」第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
誰在線上
正在瀏覽這個版面的使用者:沒有註冊會員 和 7 位訪客