天体力学/开普勒问题

天体力学

开普勒问题

开普勒问题是一个二体问题，用于找到轨道体的在一段时间内的位置和速度。给定质量、位置和速度，或者可以给出六个开普勒轨道元素。

使用开普勒轨道元素

假设轨道由半长轴 $a$ 、偏心率 $e$ 、倾角 $i$ 、升交点赤经 $\Omega$ 和近心点幅角 $\omega$ 定义。则可以找到轨道体在时间 $t$ 的位置。步骤如下

1. 计算平均近心点角

$M={\sqrt {\frac {\mu }{a^{3}}}}(t-\tau )$

其中， $\mu$ 是引力参数， $\tau$ 是轨道体经过近心点的时间。为方便起见，平均近心点角可以在 $[0,2\pi )$ 区间内表示。

2. 使用开普勒方程计算偏心近心点角

$M=E-e\sin {(E)}$

虽然开普勒方程很容易求解时间，但逆问题没有一般解。要确定给定时间点的偏心近心点角（以及航天器的位置），通常使用迭代数值方法，例如牛顿法

x_{k+1}=x_{k}-{\frac {F(x_{k})}{F'(x_{k})}}

其中

F(E)=E-e\sin(E)-M(t)

迭代采用以下形式

E_{k+1}=E_{k}-{\frac {E_{k}-e\sin(E_{k})-M(t)}{1-e\cos(E_{k})}}

对于大多数椭圆轨道，初始猜测E₀ = M 就足够了；对于偏心率大于 0.8 的轨道，可以使用E₀ = π。更好的初始猜测是可能的，但通常不需要。迭代过程重复进行，直到达到所需的精度条件，例如

|E_{k+1}-E_{k}|<\delta

其中 $\delta$ 是所需的精度。

此外，与平均异常一样，偏心异常可以限制在区间 $[0,2\pi )$ 内。

3. 计算真近点角

$\tan {\frac {\nu }{2}}={\sqrt {\frac {1+e}{1-e}}}\tan {\frac {E}{2}}$

注意，由于正切函数及其反函数的取值范围和定义域，需要进行象限检查。如果偏心异常被限制在区间 $[0,2\pi )$ 内，那么当 $E<\pi$ 时，不需要进行象限检查。

使用笛卡尔坐标

${\ddot {\boldsymbol {r}}}=-{\frac {\mu }{r^{2}}}{\hat {\boldsymbol {r}}}$

从二体问题的特例可以看出，运动的控制方程如上所示，其中 $\mu$ 是引力参数。给定一个初始位置向量， ${\boldsymbol {r}}=(x,y,z)^{\dotplus }$ ，以及对应的速度向量， ${\dot {\boldsymbol {r}}}=(u,v,w)^{\dotplus }$ ，可以使用数值积分求解未来位置。为了使用 ODE 朗格-库塔方法，控制方程必须表示为一阶 ODE 系统，如下所示。注意，下面显示的系统是矢量化控制方程的标量表示。

${\begin{cases}{\frac {dx}{dt}}=u\\{\frac {dy}{dt}}=v\\{\frac {dz}{dt}}=w\\{\frac {du}{dt}}=-{\frac {\mu x}{r^{3}}}\\{\frac {dv}{dt}}=-{\frac {\mu y}{r^{3}}}\\{\frac {dw}{dt}}=-{\frac {\mu z}{r^{3}}}\end{cases}}$

龙格-库塔方法是一种迭代半隐式数值方法，许多软件包括这种方法，例如 MATLAB 和 Python 的 SciPy 积分方法。这些方法使用自适应步长，通过使用两种方法阶数来估计局部截断误差。例如，RK45 将使用 4 阶和 5 阶，RK89 将使用 8 阶和 9 阶。误差将与用户定义的误差进行比较和检查，并根据误差是否过小或过大来调整步长，使其变小或变大。一个简单的常数正步长， $h$ ，4 阶龙格-库塔方法如下所示。

假设一个 ODE 系统如上所示，带有初始条件。

${\frac {d{\boldsymbol {y}}}{dt}}={\boldsymbol {f}}(t,{\boldsymbol {y}}),\qquad {\boldsymbol {y}}(t_{0})={\boldsymbol {y}}_{0}$

那么 4 阶龙格-库塔方法可以表述为

${\boldsymbol {y}}_{n+1}={\boldsymbol {y}}_{n}+{\frac {h}{6}}({\boldsymbol {k}}_{1}+2{\boldsymbol {k}}_{2}+2{\boldsymbol {k}}_{3}+{\boldsymbol {k}}_{4})$