二体问题的轨道模拟

二体问题的轨道模拟

为了让大家能够查询到“天体力学”方面的内容,同时锻炼我的表达和计算能力,BoJone构思了《方程与宇宙》这个主题,主要是写一些关于使用数学相对深入地讨论一些天文问题。其实我一直觉得,不用公式是无法完美地描述科学的(当然也不能纯公式),我记得霍金的《时间简史》以及《果壳中的宇宙》等之类的书,都力求不用或者尽可能少用数学公式来表达自己的观点。这种模式对于对于公众来说是很好的,但是对于希望深入研究的朋友来说却难以进行。所以我主张:宇宙是算出来的!

这个主题每一个字都是由BoJone敲击出来的,其中包括引用了《天体力学引论》里面的一些内容,以及加入了BoJone个人的一些见解。由于篇幅长及时间有限问题,BoJone打算分若干次撰写发布,并且尽可能写得通俗一点,力求让有一点微积分基础的朋友就可以弄懂。这里首先发布第一部分。由于时间匆忙等原因,可能会出现一些疏忽,欢迎大家挑错!

二体问题使假设只有两个天体,不考虑其他天体的干扰,在万有引力(牛顿经典力学)作用下如何运动的问题。这是天体力学中最简单的一类问题,也是目前唯一能够精确求解的类型。如果考虑三个或者更多天体的作用,那么只能够用到近似分析或者数值算法。

设只有太阳S和某一个行星P,建立任一空间坐标系O-xyz,$\boldsymbol{r}_s$是太阳的位置向量,$\boldsymbol{r}_p$是行星的位置向量。$\boldsymbol{r}$是行星相对于太阳的位置向量。以M和m表示太阳和行星的质量。太阳受到行星的引力
$$\boldsymbol{F}_s=G\frac{Mm}{r^2}\frac{\boldsymbol{r}}{|\boldsymbol{r}|}=G\frac{Mm}{r^3}\boldsymbol{r}$$
同时行星也受到太阳的引力
$$\boldsymbol{F}_p=G\frac{Mm}{r^2}\frac{-\boldsymbol{r}}{|\boldsymbol{r}|}=-G\frac{Mm}{r^3}\boldsymbol{r}$$
由牛顿第二定律:太阳和行星都在力的作用下产生运动
$$\boldsymbol{F}=\frac{d(mv)}{dt}=m\frac{d^2 \boldsymbol{r}}{dt^2}$$
应当有
$$\begin{aligned}\boldsymbol{F}_s=G\frac{Mm}{r^3}\boldsymbol{r}=M\frac{d^2 \boldsymbol{r}_s}{dt^2} \\ \boldsymbol{F}_p=-G\frac{Mm}{r^3}\boldsymbol{r}=m\frac{d^2 \boldsymbol{r}_p}{dt^2}\end{aligned}$$
因为$\boldsymbol{r}=\boldsymbol{r}_p-\boldsymbol{r}_s$,所以有
$$\frac{d^2 \boldsymbol{r}}{dt^2}=\frac{d^2}{dt^2}(\boldsymbol{r}_p-\boldsymbol{r}_s)=-\frac{\mu}{r^3}\boldsymbol{r}$$
其中$\mu=G(M+m)$。如果建立以太阳为原点的x-y-z直角坐标系,可以获得二体问题的微分方程组为:

$$\begin{aligned}\frac{d^2 x}{dt^2}=-\frac{\mu}{r^3}x \\ \frac{d^2 y}{dt^2}=-\frac{\mu}{r^3}y \\ \frac{d^2 z}{dt^2}=-\frac{\mu}{r^3}z\end{aligned}$$

式中$r^2=x^2+y^2+z^2$。下面我们来看一下这道方程组的解。在理论力学可以得知,行星的运动是有心力的作用,在有心力作用下的运动都是平面运动,当然我们可以从上述微分方程证实。由微分方程后两式可得:
$$y\ddot z-z\ddot y=0\Rightarrow \frac{d}{dt}(y\dot z-z\dot y)=0$$
两边积分:
$$y\dot z-z\dot y=A\tag{1}$$其中A为积分常数,同理可以得到:
$$z\dot x-x\dot z=B\tag{2}$$$$x\dot y-y\dot x=C\tag{3}$$(1),(2),(3)被称为“动量矩积分”

用x,y,z分别乘以(1),(2),(3)式后相加,得到:$Ax+By+Cz=0$
这是一个通过原点(太阳)的平面方程,这表明行星和太阳始终在同一个平面上!所以我们可以只考虑O-xy平面的方程,即
$\ddot x=-\frac{\mu x}{r^3}$,$$\ddot y=-\frac{\mu y}{r^3}\tag{4}$$动量矩积分变为一个:$$x\dot y-y\dot x=h\tag{5}$$
现在我们用极坐标讨论:设$x=rcos\theta,y=rsin\theta$,将(4)变换成
$$\ddot{r} -r\dot{\theta}^2=-\frac{\mu}{r^2}\tag{6}$$(5)变换成:
$$r^2 \dot{\theta}=h\tag{7}$$(7)式称为“面积积分”h就是向径所扫过的面积速度的两倍。由于面积速度一定,所以当然在相同时间内所扫过的面积相等啦,于是我们就证明了“开普勒第二定律”

附:这里简述一下变换过程:

设$x=r*cos\theta,y=r*sin\theta$........(00)
有$x^2+y^2=r^2$代入原方程
$\ddot{x}=-\frac{\mu cos\theta}{r^2}$,$\ddot{y}=-\frac{\mu sin\theta}{r^2}$........(0)
对(00)关于t求导
$\dot{x}=\dot{r}*cos\theta-\dot{\theta}*r*sin\theta$..........(01)
$\dot{y}=\dot{r}*sin\theta+\dot{\theta}*r*cos\theta$..........(02)
(01),(02)再次对t求导
$\ddot{x}=\ddot{r}*cos\theta-2*\dot{r}*\dot{\theta}*sin\theta-\ddot{\theta}*r*sin\theta-\dot{\theta}^2*r*cos\theta$.....................(03)
$\ddot{y}=\ddot{r}*sin\theta+2*\dot{r}*\dot{\theta}*cos\theta+\ddot{\theta}*r*cos\theta-\dot{\theta}^2*r*sin\theta.$....................(04)
$(03)*cos\theta+(04)*sin\theta$得
$\ddot{x}cos\theta+\ddot{y}sin\theta=\ddot{r}-\dot{\theta}^2*r$.....(05)
又由(0)有
$\ddot{x}cos\theta+\ddot{y}sin\theta=-{\mu}/r^2$.......(06)
对比(05),(06)即有
$$\ddot{r}-\dot{\theta}^2\cdot r=-{\mu}/r^2$$

至于从(5)式变换成(7)式,就不需要那么多技巧性的东西,直接代入求导数,然后合并就行了。

现在我们来求解由(6),(7)构成的方程组,如果希望找到行星轨道的曲线类型,那么我们需要找到r和$\theta$之间的关系式;如果我们希望计算某时刻的行星位置,那么我们得找到r或$\theta$关于时间t的关系式。

首先我们寻求轨道的曲线类型,这里我直接摘录《天体力学引论》里面的内容。
令$u=1/r$,则(7)变为:$\dot{\theta}=hu^2$,并且有
$$\begin{aligned}\dot{r}=\frac{dr}{d\theta}\frac{d\theta}{dt}=\frac{d(1/u)}{d\theta}\dot{\theta}=-1/u^2\cdot \frac{du}{d\theta}\cdot hu^2=-h\frac{du}{d\theta} \\ \ddot{r}=-h\frac{d}{d\theta}(\frac{du}{d\theta})\dot{\theta}=-h^2 u^2\frac{d^2 u}{d\theta^2}\end{aligned}$$
代入(6)得到
$$\frac{d^2 u}{d\theta^2}+u=\frac{\mu}{h^2}$$
这是一个二阶线性微分方程。它的通解是:
$$u=\frac{\mu}{h^2}[1+e \cos(\theta-\omega)]\tag{8}$$其中e和$\omega$是待定常数,还原为r,则变为:
$$r=\frac{h^2//\mu}{1+e \cos(\theta-\omega)}\tag{9}$$
由解析几何可以知道,这是一条以焦点为原点的圆锥曲线(椭圆、双曲线、抛物线,可以参考维基百科)。至此我们证明了“开普勒第一定律”,而且更加广泛(不仅包括椭圆,还包括双曲线、抛物线)。由此可以推出$h^2=\mu a(1-e^2)$,a就是圆锥曲线的半长轴,而e则是偏心率,并且$\theta=\omega$时,r最小,即行星在近日点,所以$\omega$为近日点角)

附:如何解这个微分方程?
具体求解线性微分方程可以参考维基百科

根据解法,我们可以求出$\frac{d^2 u}{d\theta^2}+u=0$的通解为$y=C_1 cos\theta+C_2 sin\theta$,加上右边的$\frac{\mu}{h^2}$即为该微分方程的通解,但是如何转换成(8)的类型呢?
不难发现其实$C_1 cos\theta+C_2 sin\theta+\frac{\mu}{h^2}$与(8)是等价的,因为
$$e \cos(\theta-\omega)=ecos\theta \cos\omega+e \sin\theta \sin\omega$$
即$ecos\omega=C_1 h^2/{\mu},esin\omega=C_2 h^2/{\mu}$

至此,我们已经完成了大部分的工作,还有最后一点点,解答二体问题就大功告成了,就是找到r或者$\theta$关于时间t的函数,才能计算某时刻的行星位置。

附:t与r的关系式
由(7)得到$\dot{\theta}=h/r^2$,代入(6)得到:
$$\ddot{r} -h^2/r^3=-\frac{\mu}{r^2}\tag{10}$$这是一个二阶微分方程,它的解很容易找出,但是这个积分太复杂:
$$\dot{r}\frac{d\dot{r}}{dr}=h^2/r^3-\frac{\mu}{r^2}$$
$\dot{r}d\dot{r}=(h^2/r^3-\frac{\mu}{r^2})dr$,两端积分
$$\dot{r}^2={2\mu}/r-h^2/r^2+K_1\tag{11}$$$$\begin{aligned}\Rightarrow {dt}/{dr}=\frac{r}{\sqrt{K_1 r^2+2\mu r-h^2}} \\ t=\int \frac{r}{\sqrt{K_1 r^2+2\mu r-h^2}}dr\end{aligned}$$
这个积分很容易的(可以参考本站的积分表),缺点是但是最终的结果形式太复杂了!

《天体力学引论》找出的是$\theta$与t的关系(即开普勒方程)......(待续)

转载到请包括本文地址:https://kexue.fm/archives/549

更详细的转载事宜请参考:《科学空间FAQ》

如果您还有什么疑惑或建议,欢迎在下方评论区继续讨论。

如果您觉得本文还不错,欢迎分享/打赏本文。打赏并非要从中获得收益,而是希望知道科学空间获得了多少读者的真心关注。当然,如果你无视它,也不会影响你的阅读。再次表示欢迎和感谢!

如果您需要引用本文,请参考:

苏剑林. (Mar. 20, 2010). 《《方程与宇宙》:二体问题的来来去去(一) 》[Blog post]. Retrieved from https://kexue.fm/archives/549

@online{kexuefm-549,
        title={《方程与宇宙》:二体问题的来来去去(一)},
        author={苏剑林},
        year={2010},
        month={Mar},
        url={\url{https://kexue.fm/archives/549}},
}