开普勒方程

<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的开普勒方程页面