這個系統模擬了兩個物體 $m_1$ 和 $m_2$ 僅受彼此間的**萬有引力**（中心力）作用下的 2D 運動。
(此處推導與彈簧力情境完全相同) ... 由於系統不受外力，總動量守恆，質心 $\vec{R}_{cm}$ 保持靜止。 $$ \vec{r_1} = \frac{m_2}{m_1 + m_2} \vec{r} \quad \text{且} \quad \vec{r_2} = -\frac{m_1}{m_1 + m_2} \vec{r} $$
這個二體問題可以被簡化為一個等效的「單體問題」，即一個質量為「折合質量」 $\mu = \frac{m_1 m_2}{m_1 + m_2}$ 的粒子，其位置向量為 $\vec{r}$，受到的力為 $\vec{F}_{12}$。 $$ \mu \ddot{\vec{r}} = \vec{F}_{12} = -G \frac{m_1 m_2}{|\vec{r}|^2} \hat{r} $$ 其中 $G$ 是萬有引力常數，$\hat{r} = \vec{r} / |\vec{r}|$ 是單位向量。
(此處推導與彈簧力情境完全相同) ... 由於萬有引力是中心力，系統的總角動量 $\vec{L}$ 守恆。 $$ L = |\vec{L}| = \mu r^2 \dot{\theta} = \text{Constant} $$ (初始 V0 的推導也相同...)
系統總能量 $E = K + U$ 也是守恆的。萬有引力的位能 $U = -G m_1 m_2 / r$。 $$ E = \frac{1}{2}\mu v^2 - \frac{G m_1 m_2}{r} $$ 將速度 $v$ 分解為徑向 $v_r = \dot{r}$ 和切向 $v_t = r \dot{\theta}$，並代入 $L = \mu r^2 \dot{\theta}$： $$ E = \underbrace{\frac{1}{2}\mu \dot{r}^2}_{\text{Radial K.E.}} + \underbrace{\left[ \frac{L^2}{2\mu r^2} - \frac{G m_1 m_2}{r} \right]}_{\text{Effective Potential } U_{eff}(r)} $$ 這就是著名的**克卜勒問題 (Kepler Problem)**。 物體的相對距離 $r$ 會在這個 $U_{eff}$ 形成的位能井中運動。根據總能量 $E$ 的不同，相對軌道 $\vec{r}(t)$ 會是圓形（$E$ 為 $U_{eff}$ 最小值）、橢圓形（$E < 0$）、拋物線形（$E = 0$）或雙曲線形（$E > 0$）。
本模擬同樣採用「尤拉-克羅默法」(Euler-Cromer Method) 進行數值積分。 我們需要計算「相對加速度 $\vec{a} = \ddot{\vec{r}}$」： $$ \vec{a} = \frac{\vec{F}_{12}}{\mu} = \left( -G \frac{m_1 m_2}{|\vec{r}|^2} \hat{r} \right) / \left( \frac{m_1 m_2}{m_1 + m_2} \right) $$ $$ \vec{a}(t) = -G \frac{(m_1 + m_2)}{|\vec{r}|^2} \hat{r}(t) $$ 在每一個微小的時間步 $\Delta t$ 中，我們更新： $$ \vec{v}(t + \Delta t) = \vec{v}(t) + \vec{a}(t) \Delta t $$ $$ \vec{r}(t + \Delta t) = \vec{r}(t) + \vec{v}(t + \Delta t) \Delta t $$ 最後再由 $\vec{r}$ 計算出 $\vec{r_1}$ 和 $\vec{r_2}$ 來繪圖。