为什么行星的运行轨道都是椭圆的?
2024-11-19 阅读 264
更新于 2024年11月21日
似乎没有计算的，我来抛砖引玉一下吧
先说结论：经典框架下，若恒星质量足够大，则行星的运行轨迹为圆锥曲线。
一、运动方程不妨设恒星位于原点，
行星运动所对应的轨迹向量为 \vec s(t)=(x(t),y(t))\vec s(t)=(x(t),y(t)) ,
求出其速度向量为 \vec{v}(t)=(x'(t),y'(t))\vec{v}(t)=(x'(t),y'(t)) ,
而加速度向量显然为 \vec{a}=(x'',y'')\vec{a}=(x'',y'')
由牛顿定理 F=ma,F=\frac{\color{blue}{GM}m}{\color{green}{r^2}}F=ma,F=\frac{\color{blue}{GM}m}{\color{green}{r^2}} 不难得到矢量等式
\vec{a}=\frac{\color{blue}k}{\color{green}{\left| s^2 \right|}}(-\frac{\vec s}{\left| s \right|})\\\vec{a}=\frac{\color{blue}k}{\color{green}{\left| s^2 \right|}}(-\frac{\vec s}{\left| s \right|})\\ 两个分量分别相等可以得到
x''=-k\frac{x}{(x^2+y^2)^{\frac{3}{2}}}\\ y''=-k\frac{y}{(x^2+y^2)^{\frac{3}{2}}}\\x''=-k\frac{x}{(x^2+y^2)^{\frac{3}{2}}}\\ y''=-k\frac{y}{(x^2+y^2)^{\frac{3}{2}}}\\ 两式相加、相除，整理一下就是
\begin{align} x''^2+y''^2&=k^2\frac{1}{(x^2+y^2)^2}&(1)\\ x''y&=xy''&(2)\\ \end{align}\\\begin{align} x''^2+y''^2&=k^2\frac{1}{(x^2+y^2)^2}&(1)\\ x''y&=xy''&(2)\\ \end{align}\\ 为了方便解这个微分方程，进行换元（这里人傻了，应该一开始就用极坐标系求解的）
x(t)=r(t)\cos(\theta(t))\\ y(t)=r(t)\sin(\theta(t))\\x(t)=r(t)\cos(\theta(t))\\ y(t)=r(t)\sin(\theta(t))\\ 对两式求二阶导
x''=r''\cos\theta-2r'\sin\theta \cdot\theta'-r\cos\theta \cdot\theta'^2-r\sin\theta\cdot\theta''\\ y''=r''\sin\theta+2r'\cos\theta \cdot\theta'-r\sin\theta \cdot\theta'^2+r\cos\theta\cdot\theta''\\x''=r''\cos\theta-2r'\sin\theta \cdot\theta'-r\cos\theta \cdot\theta'^2-r\sin\theta\cdot\theta''\\ y''=r''\sin\theta+2r'\cos\theta \cdot\theta'-r\sin\theta \cdot\theta'^2+r\cos\theta\cdot\theta''\\ 代入 (2)(2) 化简得到
2r'\theta'+r\theta''=0\\2r'\theta'+r\theta''=0\\ 利用这个式子化简微分
x''=(r''-r\theta'^2)\cos \theta\\ y''=(r''-r\theta'^2)\sin \theta\\x''=(r''-r\theta'^2)\cos \theta\\ y''=(r''-r\theta'^2)\sin \theta\\ 代入 (1)(1)
(r''-r\theta'^2)^2=\frac{k^2}{r^4}\\(r''-r\theta'^2)^2=\frac{k^2}{r^4}\\ 最后得到了一个微分方程组
\begin{cases} 2r'\theta'+r\theta''=0\\ r''-r\theta'^2=\frac{k}{r^2}\\ \end{cases}\\\begin{cases} 2r'\theta'+r\theta''=0\\ r''-r\theta'^2=\frac{k}{r^2}\\ \end{cases}\\
二、方程求解 第一个解出两者的关系
\begin{align} 2r'r\theta'+r^2\theta''&=0\\ \frac{d(r^2\theta')}{dt}&=0\\ r^2\theta'&=C \end{align}\\\begin{align} 2r'r\theta'+r^2\theta''&=0\\ \frac{d(r^2\theta')}{dt}&=0\\ r^2\theta'&=C \end{align}\\ 这个其实就是著名的开普勒第二定律，可以写作
\frac{1}{2}r^2d\theta=Cdt\\\frac{1}{2}r^2d\theta=Cdt\\
即相等时间间隔，行星与太阳的连线扫过的面积相等
solve problem NO.1,now to problem NO.2
有了关系一之后，可以将原表达式转换为不含 tt 的形式
视 rr 为 \theta\theta 的函数，并令 u=\frac{1}{r}u=\frac{1}{r}
\begin{align} r'&=\frac{dr}{d\theta}\frac{d\theta}{dt}=\frac{dr}{d\theta}\frac{C}{r^2}\\ &=-C\frac{du}{d\theta}\\ r''&=-C\frac{d^2u}{d\theta^2}\frac{d\theta}{dt}=-C^2\frac{d^2u}{d\theta^2}u^2 \end{align}\\\begin{align} r'&=\frac{dr}{d\theta}\frac{d\theta}{dt}=\frac{dr}{d\theta}\frac{C}{r^2}\\ &=-C\frac{du}{d\theta}\\ r''&=-C\frac{d^2u}{d\theta^2}\frac{d\theta}{dt}=-C^2\frac{d^2u}{d\theta^2}u^2 \end{align}\\ 代入第二个等式，
\begin{align} -C^2\frac{d^2u}{d\theta^2}u^2-C^2u^3&=ku^2\\ \frac{d^2u}{d\theta^2}+u&=A \end{align}\\\begin{align} -C^2\frac{d^2u}{d\theta^2}u^2-C^2u^3&=ku^2\\ \frac{d^2u}{d\theta^2}+u&=A \end{align}\\ 得到一个非常清爽的二阶常系数常微分方程
通解显然为
u=A+B\cos(\theta-\theta_0)\\u=A+B\cos(\theta-\theta_0)\\ 即
r=\frac{p}{1-e\cos(\theta-\theta_0)}\\r=\frac{p}{1-e\cos(\theta-\theta_0)}\\ 这是极坐标下的圆锥曲线一般方程,但是仍需要一些办法来求解未知的参数
三、初值求解根据能量守恒，引力势能和行星动能守恒，有
\begin{align} \frac{1}{2}mv_2^2-\frac{1}{2}mv_1^2=-\int_{r_1}^{r_2}\frac{GMm}{r^2}\\ v_2^2-v_1^2=2GM(\frac{1}{r_2}-\frac{1}{r_1}) \end{align}\\\begin{align} \frac{1}{2}mv_2^2-\frac{1}{2}mv_1^2=-\int_{r_1}^{r_2}\frac{GMm}{r^2}\\ v_2^2-v_1^2=2GM(\frac{1}{r_2}-\frac{1}{r_1}) \end{align}\\
另外，还有行星的曲率半径关系，
F=\frac{GMm}{r^2}=m\frac{v^2}{R}\\F=\frac{GMm}{r^2}=m\frac{v^2}{R}\\ 以下记 r'=\frac{dr}{d\theta}r'=\frac{dr}{d\theta} ,那么曲率半径就是
R=\frac{(r^2+r'^2)^{\frac{3}{2}}}{\left| r^2+2r'^2-rr'' \right|}=\frac{v^2R^2}{GM}\\R=\frac{(r^2+r'^2)^{\frac{3}{2}}}{\left| r^2+2r'^2-rr'' \right|}=\frac{v^2R^2}{GM}\\ 利用这些关系式，再加上 t_0t_0 时的初值 (r_0,\theta_0),v_0,GM(r_0,\theta_0),v_0,GM ，可以解出相关的参数
而如果行星轨道为椭圆，有
v=\sqrt{GM(\frac{2}{r}-\frac{1}{a})}\\v=\sqrt{GM(\frac{2}{r}-\frac{1}{a})}\\ 其中， aa 是半长轴，利用这个式子确定初值会方便一些。
附录1.圆锥曲线极坐标
简单添加几句关于圆锥曲线极坐标统一方程的问题
圆锥曲线的统一定义：到定点 FF 的距离与到此点外定直线 ll 的距离之比为 ee 的点的集合
其中， 01e>1 为双曲线。
在极坐标下，若 FF 位于原点，若设点 FF 到直线 ll 的距离为 pp ,则圆锥曲线的方程可以统一写作
\rho=\frac{ep}{1-e\cos \theta}\\\rho=\frac{ep}{1-e\cos \theta}\\ 2.水星进动等问题
本文主要讲的是经典物理体系下，大质量二体运动模型的解，即只可能是圆锥曲线（带退化），而行星轨道为什么都是椭圆的，主要是由于恒星质量往往非常大，且不同轨道之间的行星相互影响很小，所以可以近似认为是二体运动，而行星的轨道要求封闭，而圆锥曲线中只有椭圆是封闭轨道，因此行星的运动轨道可以认为是椭圆。
然而，很早以前，人们就观测到了水星的进动现象，
可以简单认为是水星运动分解为绕太阳的椭圆运动与椭圆轴绕太阳旋转的叠加，这得等到广义相对论出来之后才能解释了。当年爱因斯坦在发表狭义相对论后不久就曾研究过这个问题，答主曾经读过那篇大胆的论文，然而如果我没有搞错的话，他的第一篇广义相对论的论文（即解释水星进动的论文）是有问题的，不过这都是后话了，扯远了哈，总之，本篇回答主要还是集中在经典领域下的。
行星的运行轨道可以是圆形，椭圆形，抛物线形，以及双曲线形，但是太阳系演化46亿年后的今天，太阳系内所有行星的轨道都是椭圆的。
任何在引力作用下运行的行星都会沿着圆锥截面运动，如下图，这些圆锥截面曲线包括圆、椭圆、抛物线和双曲线。同时，轨道的焦点/圆点是恒星。
圆锥截面曲线这四种轨道的不同之处在于它们的偏心率。圆轨道的离心率是0，一个完美的圆轨道；椭圆轨道离心率在0和1之间，也是绝大多数行星的轨道形状；抛物线轨道离心率为1；双曲线轨道离心率大于1，这两类轨道对应的天体从太阳系外经过不会返回。
四种轨道有了这四种行星轨道的选项，我们看看为什么行星的运行轨道绝大多数都是椭圆。
首先，太阳系有 46 亿年的历史，任何具有抛物线或双曲线轨道的行星都早已不复存在了，这些轨道上的行星或者彗星，由于轨道不闭合，对于太阳系来讲都是过客；其次， 圆形轨道需要实现正好为零的偏心率，这很难，但椭圆轨道的偏心率可以在 0 到 1 之间的任何位置，这就很简单。所以，我们看到的绝大多数的行星轨道都是椭圆。
但实际上，太阳系的很多行星的轨道离心率都很小，基本上也可以近似为圆形，但严格讲还是椭圆。比如：地球离心率是0.0167，近似为圆形；离心率最小的是金星，只有0.0067，最大的是水星，为0.2056，所以水星的轨道是太阳系中最椭圆的。
牛顿发现主导行星运动的力是一个只与距离平方反比相关的定律。当时牛顿的证明是很不容易的。即使到今天，用现代微积分的做法推导出来也是不容易的，但标准的理论力学书一般都有推导。
1.我们可以先来证明，如果行星的运动轨道是椭圆，那么力是平方反比的。这里还要用到开普勒第二定律或者角动量守恒，也即
r^2\dot{\phi}=Cr^2\dot{\phi}=C （1）
另外还要用到机械能守恒，也即
\frac{1}{2}m(\dot{r}^2+r^2\dot{\phi}^2)+V(r)=E\frac{1}{2}m(\dot{r}^2+r^2\dot{\phi}^2)+V(r)=E （2）
将上式对时间求导的第二项用（1）来表示
\frac{1}{2}m(\dot{r}^2+\frac{C^2}{r^2})+V(r)=E\frac{1}{2}m(\dot{r}^2+\frac{C^2}{r^2})+V(r)=E （3）
于是我们可以发现，通过方程（1）和（3）就可以消去时间，得到轨道微分方程，于是
\frac{dr}{d\phi}=\frac{r^2}{C^2}\sqrt{E-V(r)-\frac{mC^2}{2r^2}}\frac{dr}{d\phi}=\frac{r^2}{C^2}\sqrt{E-V(r)-\frac{mC^2}{2r^2}} （4）
如果已经知道了行星运动轨道是椭圆，那么将椭圆的轨道方程代入上式中，
r=\frac{p}{1+e\cos\phi}r=\frac{p}{1+e\cos\phi} （5）
于是可以得到
V(r)\sim 1/rV(r)\sim 1/r （6）
那么力的表达式为
F=-\frac{dV}{dr}\sim-\frac{1}{r^2}F=-\frac{dV}{dr}\sim-\frac{1}{r^2} （7）
这样就证明了如果行星运动轨道是椭圆，那么所受到的力是平方反比的。细心的朋友可能发现（5）式其实不止表示椭圆，如果是抛物线和双曲线这种圆锥曲线都是可以的。因此如果轨道是抛物线和双曲线，那么所受到的力也是平方反比的。这种方法可以称之为归纳法。
2.然后再来证明，如果主导行星与太阳之间的相互作用是平方反比的，那么行星的运行轨道是椭圆。这种情况下，力是平方反比的，那么势能可以写为
V(r)=-a/r V(r)=-a/r
代入到公式（4）中，通过积分就可以导出轨道方程为（5），当然这个积分也不是直接就可以积分出来，需要先做个代换，将 rr 换成 1/\rho1/\rho ，这样的积分就很容易了，积分出来的结果是三角函数，整理下就得到（5）式了。这个演绎过程更符合题主问题的答案。这个方法称之为演绎法。
这个问题的一般解法是比内方程：
不失一般性，考虑质点在球对称中心力场 \vec F=-\frac{k}{r^2}\vec{e}_r\vec F=-\frac{k}{r^2}\vec{e}_r 中的径向运动方程。对引力问题而言， k=k_{\rm G}m_0mk=k_{\rm G}m_0m
a_{r} = \ddot{r}-r\dot{\theta}^2=-\frac{k}{r^2}a_{r} = \ddot{r}-r\dot{\theta}^2=-\frac{k}{r^2}
由于在中心力场中运动不难证明角动量 mr^2\dot{\theta}=\Thetamr^2\dot{\theta}=\Theta 是守恒量，代入上式，将 \frac{1}{r}\frac{1}{r} 作为自变量代换可以获得比内方程，然后直接求解得到椭圆轨道。这种方法楼上 @艺多不精 已经推导过了。
这里介绍一种不需要求解微分方程的办法，对高中低年级的竞赛生可能更加友好。
首先构造一个特殊的矢量 \vec{\xi}=\vec{p} \times \vec{\Theta}-m k \vec{e}_r .\vec{\xi}=\vec{p} \times \vec{\Theta}-m k \vec{e}_r . 这个矢量被称为Laplace–Runge–Lenz矢量（不要问我怎么构造出来的，问就是做梦梦到的）。我们证明他是一个守恒量。
\begin{align*} \frac{\mathrm{d} \vec{\xi}}{\mathrm{d} t} & =\frac{\mathrm{d}}{\mathrm{d} t}(\vec{p} \times \vec{\Theta})-\frac{\mathrm{d}}{\mathrm{d} t}(m k \vec{e}_r) \\ & =\frac{\mathrm{d} \vec{p}}{\mathrm{d} t} \times \vec{\Theta}+\vec{p} \times \frac{\mathrm{d} \vec{\Theta}}{\mathrm{d} t}-m k \frac{d}{\mathrm{d} t} \frac{\vec{r}}{r} \\ & =\vec{F} \times(\vec{r} \times \vec{p})+0-m k \frac{\mathrm{d}}{\mathrm{d} t} \frac{\vec{r}}{r} \\ & =-\frac{k}{r^2} \frac{\vec{r}}{r} \times(\vec{r} \times \vec{p})-m k \frac{\mathrm{d}}{\mathrm{d} t}\frac{\vec{r}}{r} \\ & =-\frac{k}{r^3}((\vec{r} \cdot \vec{p}) \vec{r}-(\vec{r} \cdot \vec{r}) \vec{p})-m k \frac{\mathrm{d}}{\mathrm{d} t}\frac{\vec{r}}{r} \\ & =-\frac{m k}{r^3}\left(\left(\vec{r} \cdot \frac{\mathrm{d} \vec{r}}{\mathrm{d} t}\right) \vec{r}-r^2 \frac{\mathrm{d} \vec{r}}{\mathrm{d} t}\right)-m k \frac{\mathrm{d}}{\mathrm{d} t}\frac{\vec{r}}{r} \\ &=-\frac{m k}{r^3}\left(r \frac{\mathrm{d} r}{\mathrm{d} t} \vec{r}-r^2 \frac{\mathrm{d} \vec{r}}{\mathrm{d} t}\right)+m k\left(\frac{1}{r^2} \frac{\mathrm{d} r}{\mathrm{d} t} \vec{r}-\frac{1}{r} \frac{\mathrm{d} \vec{r}}{\mathrm{d} t}\right)\\ & =-\frac{m k}{r^2} \frac{\mathrm{d} r}{\mathrm{d} t} \vec{r}+\frac{m k}{r} \frac{\mathrm{d} \vec{r}}{\mathrm{d} t} +\frac{m k}{r^2} \frac{\mathrm{d} r}{\mathrm{d} t} \vec{r}-\frac{m k}{r} \frac{\mathrm{d} \vec{r}}{\mathrm{d} t} \\ &=0 \end{align*}\begin{align*} \frac{\mathrm{d} \vec{\xi}}{\mathrm{d} t} & =\frac{\mathrm{d}}{\mathrm{d} t}(\vec{p} \times \vec{\Theta})-\frac{\mathrm{d}}{\mathrm{d} t}(m k \vec{e}_r) \\ & =\frac{\mathrm{d} \vec{p}}{\mathrm{d} t} \times \vec{\Theta}+\vec{p} \times \frac{\mathrm{d} \vec{\Theta}}{\mathrm{d} t}-m k \frac{d}{\mathrm{d} t} \frac{\vec{r}}{r} \\ & =\vec{F} \times(\vec{r} \times \vec{p})+0-m k \frac{\mathrm{d}}{\mathrm{d} t} \frac{\vec{r}}{r} \\ & =-\frac{k}{r^2} \frac{\vec{r}}{r} \times(\vec{r} \times \vec{p})-m k \frac{\mathrm{d}}{\mathrm{d} t}\frac{\vec{r}}{r} \\ & =-\frac{k}{r^3}((\vec{r} \cdot \vec{p}) \vec{r}-(\vec{r} \cdot \vec{r}) \vec{p})-m k \frac{\mathrm{d}}{\mathrm{d} t}\frac{\vec{r}}{r} \\ & =-\frac{m k}{r^3}\left(\left(\vec{r} \cdot \frac{\mathrm{d} \vec{r}}{\mathrm{d} t}\right) \vec{r}-r^2 \frac{\mathrm{d} \vec{r}}{\mathrm{d} t}\right)-m k \frac{\mathrm{d}}{\mathrm{d} t}\frac{\vec{r}}{r} \\ &=-\frac{m k}{r^3}\left(r \frac{\mathrm{d} r}{\mathrm{d} t} \vec{r}-r^2 \frac{\mathrm{d} \vec{r}}{\mathrm{d} t}\right)+m k\left(\frac{1}{r^2} \frac{\mathrm{d} r}{\mathrm{d} t} \vec{r}-\frac{1}{r} \frac{\mathrm{d} \vec{r}}{\mathrm{d} t}\right)\\ & =-\frac{m k}{r^2} \frac{\mathrm{d} r}{\mathrm{d} t} \vec{r}+\frac{m k}{r} \frac{\mathrm{d} \vec{r}}{\mathrm{d} t} +\frac{m k}{r^2} \frac{\mathrm{d} r}{\mathrm{d} t} \vec{r}-\frac{m k}{r} \frac{\mathrm{d} \vec{r}}{\mathrm{d} t} \\ &=0 \end{align*}
对，就是这么神奇，这个不知道怎么构造出来的量是守恒量。
于是可以根据其在某个特殊时刻的大小和方向确定这个矢量本身。比如说可以猜出来这个运动（以及其时间反演得到的运动）的轨道上有一个点距离原点最近（否则这意味着该运动与原点的距离可以任意接近，从而势能为无穷大，与运动的能量守恒矛盾），在这个点分析 \vec p\vec p 和 \vec \Theta\vec \Theta 的关系，就能得到这个矢量从坐标原点指向这个最近的一点。
现在我们向刚才那个点（近拱点）的反向建立极坐标轴。于是
\begin{align*} & \vec\xi \cdot \vec{r} = \left(\vec{p} \times \vec{\Theta}-m k \vec{e}_r\right) \cdot \vec{r}=-\xi r {\rm Cos} [\theta] \\ & \vec{r}\cdot(\vec{p} \times \vec{\Theta})-m k r=-\xi r {\rm Cos} [\theta] \\ & (\vec{r} \times \vec{p}) \cdot \vec{\Theta}-m k r=-\xi r {\rm Cos} [\theta] \text {. } \\ & \vec{\Theta} \cdot \vec{\Theta}-m k r=-\xi r {\rm Cos} [\theta] \\ & (m k-\xi \cos [\theta]) r=\Theta^2 \end{align*}\begin{align*} & \vec\xi \cdot \vec{r} = \left(\vec{p} \times \vec{\Theta}-m k \vec{e}_r\right) \cdot \vec{r}=-\xi r {\rm Cos} [\theta] \\ & \vec{r}\cdot(\vec{p} \times \vec{\Theta})-m k r=-\xi r {\rm Cos} [\theta] \\ & (\vec{r} \times \vec{p}) \cdot \vec{\Theta}-m k r=-\xi r {\rm Cos} [\theta] \text {. } \\ & \vec{\Theta} \cdot \vec{\Theta}-m k r=-\xi r {\rm Cos} [\theta] \\ & (m k-\xi \cos [\theta]) r=\Theta^2 \end{align*}
移项即得
r=\frac{\Theta}{m k-\xi \cos [\theta] }r=\frac{\Theta}{m k-\xi \cos [\theta] }
与圆锥曲线的极坐标方程 r = \frac{\beta_{\rm e}l_{\rm p}}{1-\beta_{\rm e}{\rm Cos[\theta]}}r = \frac{\beta_{\rm e}l_{\rm p}}{1-\beta_{\rm e}{\rm Cos[\theta]}} 对比可知其作圆锥曲线运动。
Q.E.D
当然可以顺便讨论一下离心率的关系，不难得到
\beta_e=\frac{\xi}{m k} \beta_e=\frac{\xi}{m k}
由于 \xi\xi 是守恒量，可以直接在近拱点求的其大小
\xi^2 =(p\Theta-mk)^2\\ =p^2 \Theta^2+m^2 k^2-2 p \Theta m k \xi^2 =(p\Theta-mk)^2\\ =p^2 \Theta^2+m^2 k^2-2 p \Theta m k
假设运动总能量为 W=\frac{p^2}{2m}-\frac{k}{r}W=\frac{p^2}{2m}-\frac{k}{r} ，由于在近拱点处速度方向垂直于位置矢径， \Theta=pr\Theta=pr . 因此，
\begin{aligned} \xi^2 & =2 m\left(W+\frac{k}{r}\right) \cdot \Theta^2+m^2 k^2-2 \frac{\Theta}{r} \Theta m k . \\ & =2 m W \Theta^2+\frac{2 m k}{r} \Theta^2+m^2 k^2-2 \frac{\Theta^2}{r} m k . \\ & =2 m W \Theta^2+m^2 k^2 . \end{aligned}\begin{aligned} \xi^2 & =2 m\left(W+\frac{k}{r}\right) \cdot \Theta^2+m^2 k^2-2 \frac{\Theta}{r} \Theta m k . \\ & =2 m W \Theta^2+\frac{2 m k}{r} \Theta^2+m^2 k^2-2 \frac{\Theta^2}{r} m k . \\ & =2 m W \Theta^2+m^2 k^2 . \end{aligned}
立即得
\begin{equation} \beta_e=\sqrt{1+\frac{2 W \Theta^2}{m k^2}} \end{equation}\begin{equation} \beta_e=\sqrt{1+\frac{2 W \Theta^2}{m k^2}} \end{equation}
