克卜勒方程式 (英語:Kepler's equation )把軌道力學 中受到連心力 影響的軌道中天體多種幾何性質聯繫起來。
圖為五種由0到1不同離心率的克卜勒方程式解
這條方程式最早由約翰內斯·克卜勒 得出,推導過程見於他在1609年出版的著作《新天文學 》的第60章[ 1] [ 2] ,而他在著作《哥白尼天文學概要 》(1621年出版)的第五卷中還提出了此方程式的一個迭代解[ 3] [ 4] 。這條方程式在物理學史和數學史上都扮演著重要的角色,特別是在古典天體力學 史上。
克卜勒方程式 如下:
M
=
E
−
e
sin
E
{\displaystyle M=E-e\sin E}
其中M 為平近點角 ,E 為偏近點角 ,而e 則為離心率 .
偏近點角E 對於計算克卜勒軌道上移動點的位置相當有用。比方說,天體在座標x = a (1 − e ) 、y = 0 和時間t = t 0 時經過近星點,那麼要找出天體在任意時間的位置的話,首先可由時間和平均運動 n 用方程式M = n (t − t 0 ) 計算出平近點角M ,然後解上述的克卜勒方程式得E ,便能由下列方程式求座標:
x
=
a
(
cos
E
−
e
)
y
=
b
sin
E
{\displaystyle {\begin{array}{lcl}x&=&a(\cos E-e)\\y&=&b\sin E\end{array}}}
其中a 為半長軸 ,而b 則為半短軸 。
由於正弦 是超越函數 ,所以克卜勒方程式是超越方程式 ,即是說不能用代數方法 解出E 。一般求E 需要用到數值分析 和級數 。
克卜勒方程式共有幾種形式。每一種形式都同一種特定軌道類型有關。標準的克卜勒方程式用於橢圓軌道(0 ≤e <1)。雙曲克卜勒方程式用於雙曲軌跡(e ≫1)。徑向克卜勒方程式用於線性(徑向)軌跡 (e =1)。拋物線軌跡(e =1)則用巴克方程式 。
當e = 0時,軌道是圓形的。增加e 會導致圓形變成橢圓。當e = 1時共有三種可能性:
拋物線軌跡;
沿著從吸引中心起始無限長射線向內或向外的軌跡;
或沿著由吸引中心和距其一定距離的點所形成的線段來回移動的軌跡。
把e 從1輕微增加會形成轉角剛好低於180度的雙曲軌道。繼續增加數值會導致轉角變小,而當e 趨向無限時,軌道則變成長度無限的直線。
雙曲克卜勒方程式如下:
M
=
e
sinh
H
−
H
{\displaystyle M=e\sinh H-H}
其中H 為雙曲偏近點角。
這條方程式是通過把M 重新定義為−1的平方根 乘上橢圓方程式的右邊而得:
M
=
i
(
E
−
e
sin
E
)
{\displaystyle M=i\left(E-e\sin E\right)}
(其中E 現為虛數),然後以iH 取代E 。
徑向克卜勒方程式如下:
t
(
x
)
=
sin
−
1
(
x
)
−
x
(
1
−
x
)
{\displaystyle t(x)=\sin ^{-1}({\sqrt {x}})-{\sqrt {x(1-x)}}}
其中t 與時間成正比,而x 則與射線上與吸引中心的距離成正比。下式是通過把克卜勒方程式成1/2 並把e 定為1而成:
t
(
x
)
=
1
2
[
E
−
sin
(
E
)
]
{\displaystyle t(x)={\frac {1}{2}}\left[E-\sin(E)\right]}
。
代入得
E
=
2
sin
−
1
(
x
)
{\displaystyle E=2\sin ^{-1}({\sqrt {x}})}
。
可以利用已知的E 值直接求出M 。但用已知M 值來求E 卻要複雜得多。因為不存在解析解 。
可以使用拉格朗日反算法 寫出克卜勒方程式解的級數 表達式,但所有e 和M 的組合都不能使級數收斂(見下文)。
文獻中對克卜勒方程式可解性的混淆已存在了四個世紀[ 5] 克卜勒本人曾對找得到通解的可能性表示懷疑。
“
對於它[克卜勒方程式]不能用演繹法求解這點我是足夠滿意的,因為弧和正弦之間有著本質上的差異。但若果我錯了的話,任何人都應該給我指出來,而那個人在我眼中就是偉大的阿波羅尼奧斯 。
”
——約翰內斯·克卜勒[ 6]
逆克卜勒方程式是所有實數值
e
{\displaystyle e}
的克卜勒方程式解:
E
=
{
∑
n
=
1
∞
M
n
3
n
!
lim
θ
→
0
+
(
d
n
−
1
d
θ
n
−
1
(
(
θ
θ
−
sin
(
θ
)
3
)
n
)
)
,
e
=
1
∑
n
=
1
∞
M
n
n
!
lim
θ
→
0
+
(
d
n
−
1
d
θ
n
−
1
(
(
θ
θ
−
e
sin
(
θ
)
)
n
)
)
,
e
≠
1
{\displaystyle E={\begin{cases}\displaystyle \sum _{n=1}^{\infty }{\frac {M^{\frac {n}{3}}}{n!}}\lim _{\theta \to 0^{+}}\!{\Bigg (}{\frac {\mathrm {d} ^{\,n-1}}{\mathrm {d} \theta ^{\,n-1}}}{\bigg (}{\bigg (}{\frac {\theta }{\sqrt[{3}]{\theta -\sin(\theta )}}}{\bigg )}^{\!\!\!n}{\bigg )}{\Bigg )},&e=1\\\displaystyle \sum _{n=1}^{\infty }{\frac {M^{n}}{n!}}\lim _{\theta \to 0^{+}}\!{\Bigg (}{\frac {\mathrm {d} ^{\,n-1}}{\mathrm {d} \theta ^{\,n-1}}}{\bigg (}{\Big (}{\frac {\theta }{\theta -e\sin(\theta )}}{\Big )}^{\!n}{\bigg )}{\Bigg )},&e\neq 1\end{cases}}}
展開得:
E
=
{
s
+
1
60
s
3
+
1
1400
s
5
+
1
25200
s
7
+
43
17248000
s
9
+
1213
7207200000
s
11
+
151439
12713500800000
s
13
+
⋯
|
s
=
(
6
M
)
1
/
3
,
e
=
1
1
1
−
e
M
−
e
(
1
−
e
)
4
M
3
3
!
+
(
9
e
2
+
e
)
(
1
−
e
)
7
M
5
5
!
−
(
225
e
3
+
54
e
2
+
e
)
(
1
−
e
)
10
M
7
7
!
+
(
11025
e
4
+
4131
e
3
+
243
e
2
+
e
)
(
1
−
e
)
13
M
9
9
!
+
⋯
,
e
≠
1
{\displaystyle E={\begin{cases}\displaystyle s+{\frac {1}{60}}s^{3}+{\frac {1}{1400}}s^{5}+{\frac {1}{25200}}s^{7}+{\frac {43}{17248000}}s^{9}+{\frac {1213}{7207200000}}s^{11}+{\frac {151439}{12713500800000}}s^{13}+\cdots {\bigg |}{s=(6M)^{1/3}},&e=1\\\\\displaystyle {\frac {1}{1-e}}M-{\frac {e}{(1-e)^{4}}}{\frac {M^{3}}{3!}}+{\frac {(9e^{2}+e)}{(1-e)^{7}}}{\frac {M^{5}}{5!}}-{\frac {(225e^{3}+54e^{2}+e)}{(1-e)^{10}}}{\frac {M^{7}}{7!}}+{\frac {(11025e^{4}+4131e^{3}+243e^{2}+e)}{(1-e)^{13}}}{\frac {M^{9}}{9!}}+\cdots ,&e\neq 1\end{cases}}}
以上的級數可用Wolfram Mathematica 的InverseSeries運算得出。
InverseSeries [ Series [ M - Sin [ M ], { M , 0 , 10 }]]
InverseSeries [ Series [ M - e Sin [ M ], { M , 0 , 10 }]]
這些函數只是簡單的麥克勞林級數 。這樣的超越函數泰勒級數寫法可被視為那些函數的定義。因此,這個解是逆克卜勒方程式的正式定義。然而,在e 非零的時候,E 並不是M 的整函數 。其導數
d
M
/
d
E
=
1
−
e
cos
E
{\displaystyle dM/dE=1-e\cos E}
在e <1時於無限複數集合中趨向零。在
E
=
±
i
cosh
−
1
(
1
/
e
)
{\displaystyle E=\pm i\cosh ^{-1}(1/e)}
有解,而在這些值時
M
=
E
−
e
sin
E
=
±
i
(
cosh
−
1
(
1
/
e
)
−
1
−
e
2
)
{\displaystyle M=E-e\sin E=\pm i\left(\cosh ^{-1}(1/e)-{\sqrt {1-e^{2}}}\right)}
(其中逆cosh取正值),而dE/dM 在這些點則趨向無限。這意味著此麥克勞林級數的收歛半徑為
cosh
−
1
(
1
/
e
)
−
1
−
e
2
{\displaystyle \cosh ^{-1}(1/e)-{\sqrt {1-e^{2}}}}
,並且當M 比這大時級數不會收歛。此級數還可用於雙曲情況,此時其收歛半徑為
cos
−
1
(
1
/
e
)
−
e
2
−
1
{\displaystyle \cos ^{-1}(1/e)-{\sqrt {e^{2}-1}}}
。當e =1時此級數只在M <2π 時收歛。
還可以寫出用e 冪表示的麥克勞林級數。這級數在e 大於拉普拉斯極限 (約為0.66)時不會收歛,與M 值無關(除非M 為2π 的倍數),但若e 小於拉普拉斯極限時則級數在任何M 值時都會收歛。級數的係數除了第一個之外(只是M 而已),都與M 成週期式關係,週期為2π 。
逆徑向克卜勒方程式(e = 1)可被寫成:
x
(
t
)
=
∑
n
=
1
∞
[
lim
r
→
0
+
(
t
2
3
n
n
!
d
n
−
1
d
r
n
−
1
(
r
n
(
3
2
(
sin
−
1
(
r
)
−
r
−
r
2
)
)
−
2
3
n
)
)
]
{\displaystyle x(t)=\sum _{n=1}^{\infty }\left[\lim _{r\to 0^{+}}\left({\frac {t^{{\frac {2}{3}}n}}{n!}}{\frac {\mathrm {d} ^{\,n-1}}{\mathrm {d} r^{\,n-1}}}\!\left(r^{n}\left({\frac {3}{2}}{\Big (}\sin ^{-1}({\sqrt {r}})-{\sqrt {r-r^{2}}}{\Big )}\right)^{\!-{\frac {2}{3}}n}\right)\right)\right]}
展開得:
x
(
t
)
=
p
−
1
5
p
2
−
3
175
p
3
−
23
7875
p
4
−
1894
3931875
p
5
−
3293
21896875
p
6
−
2418092
62077640625
p
7
−
⋯
|
p
=
(
3
2
t
)
2
/
3
{\displaystyle x(t)=p-{\frac {1}{5}}p^{2}-{\frac {3}{175}}p^{3}-{\frac {23}{7875}}p^{4}-{\frac {1894}{3931875}}p^{5}-{\frac {3293}{21896875}}p^{6}-{\frac {2418092}{62077640625}}p^{7}-\ \cdots \ {\bigg |}{p=\left({\tfrac {3}{2}}t\right)^{2/3}}}
上述結果可用Wolfram Mathematica 求得:
InverseSeries [ Series [ ArcSin [ Sqrt [ t ]] - Sqrt [( 1 - t ) t ], { t , 0 , 15 }]]
逆問題在大部份的應用中都能以數值方法求得函數的根 :
f
(
E
)
=
E
−
e
sin
(
E
)
−
M
(
t
)
{\displaystyle f(E)=E-e\sin(E)-M(t)}
可以經牛頓法 來進行迭代:
E
n
+
1
=
E
n
−
f
(
E
n
)
f
′
(
E
n
)
=
E
n
−
E
n
−
e
sin
(
E
n
)
−
M
(
t
)
1
−
e
cos
(
E
n
)
{\displaystyle E_{n+1}=E_{n}-{\frac {f(E_{n})}{f'(E_{n})}}=E_{n}-{\frac {E_{n}-e\sin(E_{n})-M(t)}{1-e\cos(E_{n})}}}
注意在這個計算E 和M 的單位是弧度角。重覆迭代直到已經達到想要的準確度(例如,當f(E) < 所需的準確度)。對大部份的橢圓軌道而言,起始值取E 0 = M (t ) 已經足夠。而對e > 0.8的軌道而言,則應取起始值{{{1}}} 。相近的手法還可以用於克卜勒方程式的雙曲形式[ 7] :66-67 。而在面對拋物線軌跡的個案時則使用巴克方程式 。
克卜勒方程式常被聲稱「沒有解釋解」;例子見這裏 (頁面存檔備份 ,存於網際網路檔案館 )。至於這是否真確取決於是否認為無限級數(或不一定收歛的級數)算解釋解。也有其他作者無理地聲稱此方程式無解;例子見M. V. K. Chari, Sheppard Joel Salon 2000 Technology & Engineering.
Danby, J. M.; Burkardt, T. M. The solution of Kepler's equation. I. Cel. Mech. 1983, 31 : 95–107. Bibcode:1983CeMec..31...95D . doi:10.1007/BF01686811 .
Conway, B. A. An improved algorithm due to Laguere for the solution of Kepler's equation. 1986. doi:10.2514/6.1986-84 .
Mikkola, Seppo. A cubic approximation for Kepler's equation. Cel. Mech. 1987, 40 (3). Bibcode:1987CeMec..40..329M . doi:10.1007/BF01235850 .
Fukushima, Toshio. A method solving kepler's equation without transcendental function evaluations. Cel. Mech. Dyn. Astron. 1996, 66 (3): 309–319. Bibcode:1996CeMDA..66..309F . doi:10.1007/BF00049384 .
Charles, E. D.; Tatum, J. B. The convergence of Newton-Raphson iteration with Kepler's equation. Cel. Mech. Dyn. Astr. 1997, 69 (4): 357–372. Bibcode:1997CeMDA..69..357C . doi:10.1023/A:1008200607490 .
Stumpf, Laura. Chaotic behaviour in the newton iterative function associated with kepler's equation. Cel. Mech. Dyn. Astr. 1999, 74 (2): 95–109. doi:10.1023/A:1008339416143 .
Palacios, M. Kepler equation and accelerated Newton method. J. Comp. Appl. Math. 2002, 138 : 335–346. Bibcode:2002JCoAM.138..335P . doi:10.1016/S0377-0427(01)00369-7 .
Boyd, John P. Rootfinding for a transcendental equation without a first guess: Polynomialization of Kepler's equation through Chebyshev polynomial equation of the sine. Appl. Num. Math. 2007, 57 (1): 12–18. doi:10.1016/j.apnum.2005.11.010 .