Answer for HW6

前三题的解析是鞠国兴解读的教材原内容。注意，它采用了一个奇怪的记号，即用 M 来表示角动量。

1 T1~T3：有心力场

1.1 前情提要

如果作用在质点上的力的方向总是沿着质点的径矢，这样的力称为有心力。如果该力的大小仅仅与质点到固定点的距离有关，则有

F = - \frac{\partial U (r)}{\partial r} e_{r} = - \frac{d U}{d r} e_{r}

在有心力场中，质点的运动位于一个固定平面内，该平面垂直于角动量 $M$ 。采用平面极坐标，有拉格朗日函数

L = \frac{m}{2} ({\dot{r}}^{2} + r^{2} {\dot{φ}}^{2}) - U (r)

$φ$ 是循环坐标，且 $L$ 不显含时间，有两个运动积分，即角动量守恒和机械能守恒

p_{φ} = M = m r^{2} \dot{φ} = const

E = \frac{1}{2} m ({\dot{r}}^{2} + r^{2} {\dot{φ}}^{2}) + U (r) = \frac{1}{2} m {\dot{r}}^{2} + \frac{1}{2} m r^{2} {(\frac{M}{m r^{2}})}^{2} + U (r) = \frac{1}{2} m {\dot{r}}^{2} + \frac{M^{2}}{2 m r^{2}} + U (r)

坐标 $r$ 与时间 $t$ 的关系为

\begin{matrix} (1) & t = \int \frac{d r}{\sqrt{\frac{2}{m} [E - U (r)] - \frac{M^{2}}{m^{2} r^{2}}}} + const \end{matrix}

运动的轨道方程为

\begin{matrix} (2) & φ = \int \frac{(M / r^{2}) d r}{\sqrt{2 m [E - U (r)] - M^{2} / r^{2}}} + const = \int \frac{(M / r^{2}) d r}{\sqrt{2 m [E - U_{eff}]}} + const \end{matrix}

其中

U_{eff} = U (r) + \frac{M^{2}}{2 m r^{2}};

称为有效势能，而 $M^{2} / (2 m r^{2})$ 称为离心势能。

$\dot{r} = 0$ 或者满足条件 $U_{eff} = E$ 的点是运动轨道的转折点，即所谓的转折点方程为

U (r) + \frac{M^{2}}{2 m r^{2}} = E .

根据转折点以及质点的运动范围, 可以分为有界运动和无界运动。有界运动通常有两个转折点。

对于 $r_{min} < r < r_{max}$ (这里 $r_{min}, r_{max}$ 是两个转折点与力心的距离) 的有界运动, 轨道可以是封闭的或开口的。轨道封闭的条件为

\begin{matrix} (3) & Δ φ = 2 \int_{r_{min}}^{r_{max}} \frac{(M / r^{2}) d r}{\sqrt{2 m (E - U (r)) - \frac{M^{2}}{r^{2}}}} \end{matrix}

等于 $2 π$ 的有理数倍, 即 $Δ φ = \frac{2 π p}{q}$ , 其中 $p$ 和 $q$ 是整数。仅对势能与 $\frac{1}{r}$ 或者 $r^{2}$ 成正比的两种有心力场, 其中有界运动的轨道是封闭的。这通常称为 Bertrand 定理。

质点可以坠落至场中心, 即 $r$ 可能趋于零的条件是

[r^{2} U (r)] |_{r \to 0} < - \frac{M^{2}}{2 m} .

1.2 T1

质点在场 $U = - α / r$ 内沿着抛物线运动，能量 $E = 0$ ，试求质点坐标对时间的依赖关系。

解：本题要求抛物线轨道相应的参数方程。由式 $(1)$ 代入 $E = 0$ ，得

t = \int \frac{r d r}{\sqrt{\frac{2 α}{m} r - \frac{M^{2}}{m^{2}}}}

作变量代换 $r = \frac{1}{2} p (1 + η^{2})$ , $p = \frac{M^{2}}{m α^{2}}$ , 则有 $d r = p η d η$ , 于是上式可变为

t = \frac{1}{2} \sqrt{\frac{m p^{3}}{α}} \int (1 + η^{2}) d η = \sqrt{\frac{m p^{3}}{α}} \frac{1}{2} η (1 + \frac{1}{3} η^{2}) .

其实也可以直接利用积分公式 $\int \frac{x d x}{\sqrt{a x - b}} = \frac{2}{3 a} (\frac{2 b}{a} + x) \sqrt{a x - b}$ 进行求解，这样就可以得到

t = \frac{m}{3 α} (\frac{M^{2}}{m α} + r) \sqrt{\frac{2 α}{m} r - \frac{M^{2}}{m^{2}}} .

作前述的变量变换可得与前面相同的结果。因为 $x = r \cos φ$ ，利用轨道方程 $p / r = 1 + \cos φ$ ，则有

x = r \cos φ = p - r = p - \frac{1}{2} p (1 + η^{2}) = \frac{1}{2} p (1 - η^{2}) .

利用关系 $x^{2} + y^{2} = r^{2}$ ，有

y = \sqrt{r^{2} - x^{2}} = \sqrt{{[\frac{1}{2} p (1 + η^{2})]}^{2} - {[\frac{1}{2} p (1 - η^{2})]}^{2}} = p η .

注意，因为 $r$ 的取值范围是 $[0, \infty)$ ，则有 $- \infty < η < \infty$ .

1.3 T2

质点在有心力场 $U = - α / r^{2} (α > 0)$ 内运动, 试求积分形式运动方程。

解：将函数代入式 $(1)$ , 得

\begin{aligned} t & = \int \frac{d r}{\sqrt{\frac{2}{m} (E - U - \frac{M^{2}}{2 m r^{2}})}} = \sqrt{\frac{m}{2}} \int \frac{d r}{\sqrt{E + \frac{α}{r^{2}} - \frac{M^{2}}{2 m r^{2}}}} \\ = \sqrt{\frac{m}{2}} \int \frac{r d r}{\sqrt{E r^{2} + (α - \frac{M^{2}}{2 m})}} = \sqrt{\frac{m}{2}} \frac{1}{E} \sqrt{E r^{2} + α - \frac{M^{2}}{2 m}} . \end{aligned}

这给出了 $r$ 与时间 $t$ 的关系 $r = r (t)$ ,

r = \sqrt{\frac{2 E}{m} t^{2} - \frac{α}{E} + \frac{M^{2}}{2 m E}}

为求 $φ$ 与时间 $t$ 的关系 $φ = φ (t)$ , 先求 $r = r (φ)$ 。由式 $(2)$ , 适当选取 $φ$ 的起始线可以使得其中的 $const = 0$ , 则有

φ = \int \frac{M d r / r^{2}}{\sqrt{2 m (E + \frac{α}{r^{2}} - \frac{M^{2}}{2 m r^{2}})}} = \int \frac{- M d u}{\sqrt{2 m [E + (α - \frac{M^{2}}{2 m}) u^{2}]}}

积分的结果取决于参数的取不同值, 可以分为几种情况.

a) 当 $E > 0$ , $\frac{M^{2}}{2 m} - α = β > 0$ 时, 则上式变为

φ = - \frac{M}{\sqrt{2 m β}} \int \frac{d u}{\sqrt{\frac{E}{β} - u^{2}}} = \frac{M}{\sqrt{2 m β}} \arccos \frac{u}{\sqrt{E / β}},

即

\frac{1}{r} = u = \sqrt{\frac{E}{β}} \cos (\frac{\sqrt{2 m β}}{M} φ) = \sqrt{\frac{2 m E}{M^{2} - 2 m α}} \cos (φ \sqrt{1 - \frac{2 m α}{M^{2}}}) .

当 $φ = 0$ 时, $r$ 有极小值

\begin{matrix} (i) & r_{min} = \sqrt{\frac{M^{2} - 2 m α}{2 m E}}, \end{matrix}

即 $φ = 0$ 所在位置是近心点与力心的连线。因为 $r > 0$ , 则 $φ$ 的取值受到限制. 当 $\cos (φ \sqrt{1 - \frac{2 m α}{M^{2}}}) = 0$ , 即 $φ = \pm \frac{π}{2} / \sqrt{1 - \frac{2 m α}{M^{2}}}$ 时, $r \to \infty$ , 故 $φ$ 的取值范围为

- \frac{π}{2} / \sqrt{1 - \frac{2 m α}{M^{2}}} < φ < \frac{π}{2} / \sqrt{1 - \frac{2 m α}{M^{2}}} .

这种情况下的典型轨道是螺线, 常称为 Cotes 螺线。图例：

r = \sec(0.05 \theta + 1)|0<\theta< 10

b) 当 $E > 0$ , $α - \frac{M^{2}}{2 m} = β > 0$ 时, 则有

φ = - \frac{M}{\sqrt{2 m β}} \int \frac{d u}{\sqrt{\frac{E}{β} + u^{2}}} = - \frac{M}{\sqrt{2 m β}} arsinh \frac{u}{\sqrt{E / β}},

即

\frac{1}{r} = u = - \sqrt{\frac{E}{β}} \sinh (\frac{\sqrt{2 m β}}{M} φ) = - \sqrt{\frac{2 m E}{2 m α - M^{2}}} \sinh (φ \sqrt{\frac{2 m α}{M^{2}} - 1}) .

同样 $φ$ 的取值也受到限制以保证 $r > 0$ , 这里 $φ < 0$ 即可.

c) 当 $E < 0$ , $α - \frac{M^{2}}{2 m} = β > 0$ 时, 则有

φ = - \frac{M}{\sqrt{2 m β}} \int \frac{d u}{\sqrt{u^{2} - \frac{| E |}{β}}} = - \frac{M}{\sqrt{2 m β}} arccosh \frac{u}{\sqrt{| E | / β}},

即有

\frac{1}{r} = u = \sqrt{\frac{2 m | E |}{2 m α - M^{2}}} \cosh (φ \sqrt{\frac{2 m α}{M^{2}} - 1}) .

当 $φ = 0$ 时, $r$ 有极大值

r_{max} = \sqrt{\frac{2 m α - M^{2}}{2 m | E |}},

即 $φ = 0$ 所在位置是远心点与力心的连线。

需要说明的是:

(i) $E < 0$ , $\frac{M^{2}}{2 m} - α = β > 0$ 这种情况需要排除, 因为这将导致复数的 $φ$ .
(ii) 在 b) 和 c) 两种情况下, 当 $| φ | \to \infty$ 时, 均有 $r \to 0$ , 即质点可以“坠落”到力心.
(iii) 本问题用 Binet 公式可以比较方便的求解, 参见本节 T3.

本问题中质点的运动特征也可以通过有效势能进行分析. 有效势能为

U_{eff} = - \frac{α}{r} + \frac{M^{2}}{2 m r^{2}} = - (α - \frac{M^{2}}{2 m} \frac{1}{r}) \frac{1}{r} .

当 $α - \frac{M^{2}}{2 m} = β > 0$ 时, 有效势能为负, 表现为吸引力. 因此, 质点将冲向力心运动并最终“坠落”到力心; 当 $α - \frac{M^{2}}{2 m} < 0$ 时, 有效势能为正, 表现为相斥力, 此时有一个转折点, 可由 $E = U_{eff} = - (α - \frac{M^{2}}{2 m}) \frac{1}{r_{min}^{2}}$ 确定, 即有

r_{min} = \sqrt{\frac{\frac{M^{2}}{2 m} - α}{E}} = \sqrt{\frac{M^{2} - 2 m α}{2 m E}} .

这与式 $(i)$ 一致, 而且显然要求 $E > 0$ .

1.4 T3

在势能 $U = - α / r$ 上增加一个小的修正 $δ U$ , 有限运动的轨道不再封闭, 并且每运动一圈轨道的近心点都有很小的角度改变量 $δ φ$ 。在下面的情况下求 $δ φ$ ： a) $δ U = β / r^{2}$ , b) $δ U = γ / r^{3}$ .

解：根据公式 $(4)$ , 当 $r$ 从 $r_{min}$ 变到 $r_{max}$ 再变到 $r_{min}$ 时, 角变量 $φ$ 的变化量 $Δ φ$ , 也即相邻近心点对力心的夹角为

\begin{matrix} (ii) & \begin{aligned} Δ φ & = 2 \int_{r_{min}}^{r_{max}} \frac{(M / r^{2}) d r}{\sqrt{2 m [E - U (r)] - M^{2} / r^{2}}} \\ = - 2 \frac{\partial}{\partial M} \int_{r_{min}}^{r_{max}} \sqrt{2 m [E - U (r)] - \frac{M^{2}}{r^{2}}} d r . \end{aligned} \end{matrix}

将 $U = - α / r + δ U$ 代入上式, 注意到 $δ U$ 为小量, 可以对上式中的被积函数作 Taylor 展开, 并保留到 $δ U$ 的一阶项, 有

\begin{aligned} \sqrt{2 m [E - U (r)] - \frac{M^{2}}{r^{2}}} & = \sqrt{2 m (E + \frac{α}{r} - δ U) - \frac{M^{2}}{r^{2}}} \\ \approx \sqrt{2 m (E + \frac{α}{r}) - \frac{M^{2}}{r^{2}}} - \frac{m δ U}{\sqrt{2 m (E + \frac{α}{r}) - \frac{M^{2}}{r^{2}}}} . \end{aligned}

由此式 $(ii)$ 变为

Δ φ = - 2 \frac{\partial}{\partial M} \int_{r_{min}}^{r_{max}} \sqrt{2 m (E + \frac{α}{r}) - \frac{M^{2}}{r^{2}}} d r + \frac{\partial}{\partial M} \int_{r_{min}}^{r_{max}} \frac{2 m δ U d r}{\sqrt{2 m (E + \frac{α}{r}) - \frac{M^{2}}{r^{2}}}} .

其中第一项对应于在势能为 $- α / r$ 的场中质点运动轨道的相邻近心点对力心的夹角, 它等于 $2 π$ 。上式中的第二项是势场修正引起的角变量的变化, 记为 $δ φ$ , 即

\begin{matrix} (iii) & δ φ = \frac{\partial}{\partial M} \int_{r_{min}}^{r_{max}} \frac{2 m δ U d r}{\sqrt{2 m (E + \frac{α}{r}) - \frac{M^{2}}{r^{2}}}} . \end{matrix}

利用力场未受修正时 $r$ 与 $φ$ 之间的关系 (参见式 $(2)$ )

d r = \frac{r^{2}}{M} \sqrt{2 m (E - U - \frac{M^{2}}{2 m r^{2}})} d φ = \frac{r^{2}}{M} \sqrt{2 m (E + \frac{α}{r} - \frac{M^{2}}{2 m r^{2}})} d φ,

则式 $(iii)$ 可以改写为

\begin{matrix} (iv) & δ φ = \frac{\partial}{\partial M} (\frac{2 m}{M} \int_{0}^{π} r^{2} δ U d φ) \end{matrix}

a) 当 $δ U = β / r^{2}$ 时, 则有

\int_{0}^{π} r^{2} δ U d φ = \int_{0}^{π} r^{2} \frac{β}{r^{2}} d φ = β \int_{0}^{π} d φ = π β,

于是由式 $(iv)$ 可得

\begin{matrix} (v) & δ φ = \frac{\partial}{\partial M} (\frac{2 m}{M} π β) = - \frac{2 m π β}{M^{2}} = - \frac{2 π β}{p α} . \end{matrix}

对于这种情况, 势能实际上为

U = - \frac{α}{r} + \frac{β}{r^{2}};

另外，质点在该力场中的运动是可以解析求解的。为此我们采用 Binet 公式 $\frac{d^{2} u (φ)}{d φ^{2}} + u = - \frac{m}{M^{2}} \frac{d U (u)}{d u}$ 求解, 用 $u = 1 / r$ 表示时, 有

U = - α u + β u^{2} .

\frac{d^{2} u}{d φ^{2}} + u = - \frac{m}{M^{2}} (- α + 2 β u),

即

\frac{d^{2} u}{d φ^{2}} + (1 + \frac{2 m β}{M^{2}}) u = \frac{m α}{M^{2}} .

这个方程的解为齐次方程的通解和非齐次方程的特解之和, 即

u = A \cos (φ \sqrt{1 + \frac{2 m β}{M^{2}}}) + \frac{m α}{M^{2} + 2 m β},

其中 $A$ 为积分常数, 这里选取合适的初始条件使得初位相为零, 由上式, 可得

\begin{matrix} (vi) & r = \frac{1}{A \cos (φ \sqrt{1 + \frac{2 m β}{M^{2}}}) + \frac{m α}{M^{2} + 2 m β}} = \frac{p^{'}}{1 + e^{'} \cos (φ \sqrt{1 + \frac{2 m β}{M^{2}}})}, \end{matrix}

其中

p^{'} = \frac{M^{2} + 2 m β}{m α}, e^{'} = \frac{A (M^{2} + 2 m β)}{m α} .

由式 $(vi)$ 可知, 轨道是有限的, 但是, 当 $φ$ 从 $0$ 变到 $2 π$ 时, $r$ 从其极小值 $r_{min} = \frac{p^{'}}{1 + e^{'}}$ 变到 $r |_{φ = 2 π} \neq r |_{φ = 0} = r_{min}$ ; 不回到初始的位置, 即轨道不封闭.

从式 $(vi)$ 可知, 相邻两个近心点 (或远心点) 之间的夹角 $Δ φ$ 满足关系

Δ φ \sqrt{1 + \frac{2 m β}{M^{2}}} = 2 π,

即

Δ φ = \frac{2 π}{\sqrt{1 + \frac{2 m β}{M^{2}}}} .

当对力场没有修正，即 $β = 0$ 时，这个夹角应为 $2 π$ 。可见为力势的变化导致近心角的改变为

δ φ = Δ φ - 2 π = 2 π [\frac{1}{\sqrt{1 + \frac{2 m β}{M^{2}}}} - 1] \approx - \frac{2 π m β}{M^{2}} .

其中最后一个等号是假定 $β$ 为小量的结果, 这与式 $(v)$ 是相同的。

b) 当 $δ U = γ / r^{3}$ 时, 有

\begin{aligned} \int_{0}^{π} r^{2} δ U d φ & = \int_{0}^{π} r^{2} \frac{γ}{r^{3}} d φ = γ \int_{0}^{π} \frac{1}{r} d φ \\ = \frac{γ}{p} \int_{0}^{π} (1 + e \cos φ) d φ \\ = \frac{γ π}{p} = \frac{γ π m α}{M^{2}}, \end{aligned}

其中 $p = \frac{M^{2}}{m α}, e = \sqrt{1 + \frac{2 E M^{2}}{m α^{2}}}$ , 将其代入式 $(iv)$ 可得

δ φ = \frac{\partial}{\partial M} (\frac{2 m}{M} \frac{γ π m α}{M^{2}}) = - \frac{6 m^{2} α γ π}{M^{4}} = - \frac{6 γ π}{p^{2} α} .

在这种情况下, 势能为

U = - \frac{α}{r} + \frac{γ}{r^{3}},

相应地, 有心力为

F = - \frac{d U}{d r} = - \frac{α}{r^{2}} + \frac{3 γ}{r^{4}},

其中第二项可以视为广义相对论对万有引力 (第一项) 的修正, 由此可以解释水星近日点的进动问题。在经典力学中该力场相关轨道的求解以及运动性质的讨论可参考有关文献。

2 T4

记

\begin{matrix} (1) & \vec{Q} = \vec{L} \times \vec{A} = \vec{L} \times (\vec{p} \times \vec{L} - μ k \hat{r}) = c o n s t \end{matrix}

其中

\vec{L} \times (\vec{p} \times \vec{L}) = \vec{p} (\vec{L} \cdot \vec{L}) - \vec{L} (\vec{L} \cdot \vec{p}) = \vec{p} L^{2} = L^{2} μ \vec{v}

代入(1)式可得到：

\vec{v} = \frac{\vec{Q}}{m L^{2}} + \frac{k}{L^{2}} (\vec{L} \times \hat{r})

右侧第一项是常矢量，且由于 $\vec{L} ⊥ \hat{r}$ ，右侧第二项的模长恒为 $\frac{k}{L}$ 。显然，第一项就是原点偏离圆心的大小，第二项就是速度图的半径大小。

3 T5：选做题

3.1 (a)

选取坐标参数如下：

$m_{3}$ 在静止系中的速度是

\vec{v} = \vec{v^{'}} + \vec{ω} \times \vec{r} = (\dot{x}, \dot{y}) + \vec{ω} \times (x, y)

于是拉氏量是

\begin{aligned} L & = \frac{1}{2} m_{3} | \vec{v^{'}} + \vec{ω} \times \vec{r} |^{2} - V, & V = - \frac{G m_{1} m_{3}}{r_{13}} - \frac{G m_{2} m_{3}}{r_{23}} \\ = \frac{1}{2} m_{3} v^{' 2} + m_{3} \vec{v^{'}} \cdot (\vec{ω} \times \vec{r}) - V_{e f f}, & V_{e f f} = V - \frac{1}{2} m_{3} | \vec{ω} \times \vec{r} |^{2} \\ = \frac{1}{2} m_{3} ({\dot{x}}^{2} + {\dot{y}}^{2}) + m_{3} ω (x \dot{y} - y \dot{x}) \end{aligned}

代入拉格朗日方程得到 $m_{3}$ 的运动方程：

\begin{array}{r} m_{3} \ddot{x} - 2 m_{3} ω \dot{y} = - \frac{\partial V_{e f f}}{\partial x} \\ m_{3} \ddot{y} + 2 m_{3} ω \dot{x} = - \frac{\partial V_{e f f}}{\partial y} \end{array}

显然上式第二项就是科里奥利力。对于平衡点，要求 $\dot{\vec{r}} = \ddot{\vec{r}} = 0$ ，即上式左侧为 0，即

\begin{aligned} (1) & 0 & = \frac{1}{m_{3}} \frac{\partial v_{eff}}{\partial x} = G m_{1} \frac{x + a μ / m_{1}}{r_{13}^{3}} + G m_{2} \frac{x - a μ / m_{2}}{r_{23}^{3}} - ω^{2} x \\ (2) & 0 & = \frac{1}{m_{3}} \frac{\partial V_{eff}}{\partial y} = G m_{1} \frac{y}{r_{13}^{3}} + G m_{2} \frac{y}{r_{23}^{3}} - ω^{2} y . \end{aligned}

i) 若 $y \neq 0$ ，(2)就意味着：

G \frac{m_{1}}{r_{13}^{3}} + G \frac{m_{2}}{r_{23}^{3}} = ω^{2}

代入(1)消去 $r_{13}$ 或 $r_{23}$ 就得到：

ω^{2} = \frac{G (m_{1} + m_{2})}{r_{13}^{3}} = \frac{G (m_{1} + m_{2})}{r_{23}^{3}}

与 $ω^{2} = \frac{G (m_{1} + m_{2})}{a^{3}}$ 对比可知

r_{13} = r_{23} = a

也就是 $m_{1}, m_{2}, m_{3}$ 构成等边三角形，对应的解就是下图的 $L_{4}$ 和 $L_{5}$ 。

ii) 若 $y = 0$ ，(1)式就简化为

G m_{1} \frac{x + a μ / m_{1}}{| x - a μ / m_{1} |^{3}} + G m_{2} \frac{x - a μ / m_{2}}{| x - a μ / m_{2} |^{3}} - ω^{2} x = 0

该方程没有解析解，只能用数值方法去逼近，它有三个解，对应于上图的 $L_{1}, L_{2}, L_{3}$ 。当然，在限制性三体问题中 $(m_{3} ≪ m_{2} ≪ m_{1})$ ，上式有近似解，就比如常见的日-地-月三者，就构成一个这样的限制性三体系统。

3.2 (b)

这道题没有给定参数的具体值，所以同学们可以任取(比如 $m_{1} = m_{2} = m_{3} = 1$ )。这里就以地-月-卫星系统为例，因为它是限制性三体系统，所以

{\begin{aligned} x_{L_{1}} = R - R {(\frac{m_{2}}{m_{2} + 3 m_{1}})}^{1 / 3} \approx 0.84 R \approx 3.2 \times 10^{5} k m \\ x_{L_{2}} = R + R {(\frac{m_{2}}{m_{2} + 3 m_{1}})}^{1 / 3} \approx 1.16 R \approx 4.5 \times 10^{5} k m \\ x_{L_{3}} = - R + R (\frac{7 m_{2}}{5 m_{2} + 12 m_{1}}) \approx - 0.99 R \approx - 3.85 \times 10^{5} k m \\ {\vec{r}}_{L_{4}}, {\vec{r}}_{L_{5}} = (\frac{1}{2} R, \pm \frac{\sqrt{3}}{2} R) \approx (1.96, \pm 3.39) \times 10^{5} k m \end{aligned}

对应的绘图如下:

至于题目中的等势面，描述的其实就是相对于旋转系的速度为 0 的零速度面，我们一般称之为 Hill 曲面。它对应的表达式是(这里还是采用和第一问一样的，以质心为原点，而不是题目中以 $m_{1}$ 为原点)：

\frac{1}{m_{3}} V_{e f f} = - \frac{1}{2} ω^{2} (x^{2} + y^{2}) - \frac{G m_{1}}{r_{13}} - \frac{G m_{2}}{r_{23}} \equiv C

绘制这个图对 C 的精度要求比较高，画出来是

Appendix 代码 %

上面两幅图用 python 绘制，代码是

# 绘制拉格朗日点
import numpy as np
import matplotlib.pyplot as plt

# 相对质量
earth_mass = 1.0
moon_mass = 7.342e22/(5.972e24)
R = 1.0   # 地月距离(归一化)

earth_position = np.array([0.0, 0.0])
moon_position = np.array([R, 0.0])
lagrange_points = [
    {"name": "L1", "position": np.array([0.84, 0.0])},
    {"name": "L2", "position": np.array([1.16, 0.0])},
    {"name": "L3", "position": np.array([-1.0, 0.0])},
    {"name": "L4", "position": np.array([0.5, np.sqrt(3)/2])},
    {"name": "L5", "position": np.array([0.5, -np.sqrt(3)/2])},
]

plt.figure(figsize=(8, 8))
plt.scatter(earth_position[0], earth_position[1], color='blue', s=100, label='Earth')
plt.scatter(moon_position[0], moon_position[1], color='gray', s=30, label='Moon')

for point in lagrange_points:
    plt.scatter(point["position"][0], point["position"][1], color='red', marker='x', s=50)
    plt.text(point["position"][0], point["position"][1], point["name"], fontsize=12, ha='right')
# 绘制月球轨道
theta = np.linspace(0, 2*np.pi, 100)
r = R * np.ones_like(theta)
plt.plot(r * np.cos(theta), r * np.sin(theta), color='gray', linestyle='-', label='Moon Orbit') # ecc=0.055忽略

plt.title('Lagrange Points', fontsize=16)
plt.xlabel('X (Normalized)', fontsize=12)
plt.ylabel('Y (Normalized)', fontsize=12)
plt.legend()
plt.grid(True)
plt.axhline(0, color='black', linewidth=0.5, linestyle='--')
plt.axvline(0, color='black', linewidth=0.5, linestyle='--')
plt.xlim(-1.5, 1.5)
plt.ylim(-1.5, 1.5)
plt.show()

# 绘制希尔曲线族
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

G = 6.6743e-11  # m^3 kg^-1 s^-2
M = 5.972e24     # kg
m = 7.342e22     # kg
mu_1 = 398600   # km^3/s^2
mu_2 = 4903.02  # km^3/s^2
R = 3.844e5      # km
omega = np.sqrt((mu_1 + mu_2) / (R**3))  # 月球绕地球的角速度

# 网格
x = np.linspace(-6e5, 6e5, 20000)
y = np.linspace(-6e5, 6e5, 20000)
z = np.array([0])
x, y, z = np.meshgrid(x, y, z, indexing='ij')

earth_x = (m / (M + m)) * R # 假设地月质心在原点
moon_x = -(M / (M + m)) * R

# spacecraft 到地球和月球的距离
r_e = np.sqrt((x - earth_x)**2 + (y)**2 + (z)**2)
r_m = np.sqrt((x - moon_x)**2 + (y)**2 + (z)**2)
C_jacobi = 2 * mu_1 / r_e + 2 * mu_2 / r_m + omega**2 * (x**2 + y**2)  # 计算势能

def plot(C, ax):
    # 绘制希尔曲面
    mask = (C_jacobi < C) & (C_jacobi > C-0.001)
    print(f"num of points: {np.sum(mask)}")
    ax.scatter(x[mask], y[mask], marker='.', s=10, alpha=1, label=f"C={C}")

    # 设置坐标轴标签等
    ax.set_xlabel('X(km)')
    ax.set_ylabel('Y(km)')
    ax.ticklabel_format(style='sci', scilimits=(0,0), axis='both')

fig, ax = plt.subplots(figsize=(12, 12))
for i, C in enumerate([3.500, 3.343, 3.338, 3.326, 3.300, 3.159, 3.150, 3.135]):
    print(f"C={C}")
    plot(C, ax=ax)
    
plt.legend()
# plt.tight_layout()
plt.show()

1 T1~T3：有心力场

1.1 前情提要

1.2 T1

1.3 T2

1.4 T3

2 T4

3 T5： 选做题

3.1 (a)

3.2 (b)

Appendix 代码 %

3 T5：选做题