克卜勒方程

<image>

克卜勒方程（Kepler's equation）把軌道力學中受到連心力影響的軌道中天體多種幾何性質聯繫起來。

這條方程最早由約翰內斯·克卜勒得出，推導過程見於他在1609年出版的著作《新天文學》的第60章，而他在著作《哥白尼天文學概要》（1621年出版）的第五卷中還提出了此方程的一個迭代解。這條方程在物理學史和數學史上都扮演著重要的角色，特別是在古典天體力學史上。

## 方程

克卜勒方程如下：

$$M = E - e \sin E$$

其中M為平近點角，E為偏近點角，而e則為離心率.

偏近點角E對於計算克卜勒軌道上移動點的位置相當有用。比方說，天體在座標x = a(1 − e)、y = 0和時間t = t0時經過近星點，那麼要找出天體在任意時間的位置的話，首先可由時間和平均運動n用方程M = n(t − t0)計算出平近點角M，然後解上述的克卜勒方程得E，便能由下列方程求座標：

$$\begin{array}{lcl} x & = & a (\cos E - e) \\ y & = & b \sin E \end{array} <pre>$$

其中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$$

其中H為雙曲偏近點角。

這條方程是通過把M重新定義為−1的平方根乘上橢圓方程的右邊而得：

$$M = i \left( E - e \sin E \right)$$

（其中E現為虛數），然後以iH取代E。

### 徑向克卜勒方程

徑向克卜勒方程如下：

$$t(x) = \sin^{-1}( \sqrt{ x } ) - \sqrt{ x ( 1 - x ) }$$

其中t與時間成正比，而x則與射線上與吸引中心的距離成正比。下式是通過把克卜勒方程成1/2並把e定為1而成：

$t(x) = \frac{1}{2}\left[ E - \sin( E ) \right]$。

代入得

$E = 2 \sin^{-1}(\sqrt{ x })$。

## 逆問題

可以利用已知的E值直接求出M。但用已知M值來求E卻要複雜得多。因為不存在解析解。

可以使用拉格朗日反算定理寫出克卜勒方程解的級數表達式，但所有e和M的組合都不能使級數收斂（見下文）。

文獻中對克卜勒方程可解性的混淆已存在了四個世紀克卜勒本人曾對找得到通解的可能性表示懷疑。

### 逆克卜勒方程

逆克卜勒方程是所有實數值$e$的克卜勒方程解：

$$<pre> E = </pre> \begin{cases} \displaystyle \sum_{n=1}^\infty <pre> {\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) </pre> \Bigg) , & e = 1 \\ \displaystyle \sum_{n=1}^\infty { \frac{ M^n }{ n! } } \lim_{\theta \to 0^+} \! \Bigg( <pre> \frac{\mathrm{d}^{\,n-1}}{\mathrm{d}\theta^{\,n-1}} \bigg( \Big( \frac{ \theta }{ \theta - e \sin(\theta)} \Big)^{\!n} \bigg) </pre> \Bigg) , & e \ne 1 \end{cases}$$

展開得：

$$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} + <pre> \frac{151439}{12713500800000 }s^{13}+ \cdots \bigg|{s = ( 6 M )^{1/3}} , & e = 1\\ </pre> \\ \displaystyle <pre> \frac{1}{1-e} M </pre> - \frac{e}{(1-e)^4 } \frac{M^3}{3!} + \frac{(9 e^2 + e)}{(1-e)^7 } \frac{M^5}{5!} - \frac{(225 e^3 + 54 e^2 + e) }{(1-e)^{10} } \frac{M^7}{7!} + \frac{ (11025 e^4 + 4131 e^3 + 243 e^2 + e) }{(1-e)^{13} } \frac{M^9}{9!}+ \cdots , & e \ne 1 \end{cases}$$

以上的級數可用Wolfram Mathematica的InverseSeries運算得出。

```mathematica
InverseSeries[Series[M - Sin[M], {M, 0, 10}]]
```

```mathematica
InverseSeries[Series[M - e Sin[M], {M, 0, 10}]]
```

這些函數只是簡單的麥克勞林級數。這樣的超越函數泰勒級數寫法可被視為那些函數的定義。因此，這個解是逆克卜勒方程的正式定義。然而，在e非零的時候，E並不是M的整函數。其導數

$$dM/dE=1-e\cos E$$

在e<1時於無限複數集合中趨向零。在$E=\pm i\cosh^{-1}(1/e)$有解，而在這些值時

$$M=E-e\sin E=\pm i\left(\cosh^{-1}(1/e)-\sqrt{1-e^2}\right)$$

（其中逆cosh取正值），而dE/dM在這些點則趨向無限。這意味著此麥克勞林級數的收歛半徑為$\cosh^{-1}(1/e)-\sqrt{1-e^2}$，並且當M比這大時級數不會收歛。此級數還可用於雙曲情況，此時其收歛半徑為$\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）可被寫成：

$$<pre> x( t ) = \sum_{n=1}^{ \infty } </pre> \left[ \lim_{ r \to 0^+ } \left( <pre> {\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 - \frac{1}{5} p^2 - \frac{3}{175}p^3 <pre> - \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求得:

```mathematica
InverseSeries[Series[ArcSin[Sqrt[t]] - Sqrt[(1 - t) t], {t, 0, 15}]]
```

## 逆問題的數值近似

逆問題在大部份的應用中都能以數值方法求得函數的根：

$$f(E) = E - e \sin(E) - M(t)$$

可以經牛頓法來進行迭代：

$$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) 0.8的軌道而言，則應取起始值E_0=π。相近的手法還可以用於克卜勒方程的雙曲形式。而在面對拋物線軌跡的個案時則使用拋物線軌跡。

### 定點迭代

另一種有關的解法由考慮$E=M+e\sin{E}$開始。重覆地將E的值代入右方式子就能得到簡單的求$E(e,M)$的定點迭代演算法。此方法與克卜勒1621年的解相同。

```mathematica
function E(e,M,n)
    E = M
    for k = 1 to n
        E = M + e*sin E
    next k
    return E
```

迭代次數$n$取決於$e$的數值。雙曲形式則用相近的$H= e \sinh H - M$。

這個方法與上文牛頓法解的關係如下：

$$E_{n+1} = E_{n} - \frac{ E_{n} - e \sin(E_{n}) - M(t) }{ 1 - e \cos(E_{n})} = E_{n} + \frac {(M + e \sin{E_{n}} - E_{n})(1 + e \cos{E_{n}})}{1 - e^2 (\cos{E_{n}})^2}$$

若$M-E_{n}$及$e$的數值小的話，取一次項得近似：

$E_{n+1} \approx M + e \sin{E_{n}}$。

## 外部連結

- Wolfram Mathworld的克卜勒方程頁面