天體力學/開普勒問題

天體力學

開普勒問題

開普勒問題是一個二體問題，用於找到軌道體的在一段時間內的位置和速度。給定質量、位置和速度，或者可以給出六個開普勒軌道元素。

使用開普勒軌道元素

假設軌道由半長軸 $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})$