图

天体运动的简单数值计算

         

预备知识 万有引力, 弹簧振子受迫运动的简单数值计算

   下面我们来用一种极其简单的算法对单个天体在中心天体的万有引力作用下的运动进行数值计算. 事实上该问题存在解析解(见开普勒三定律), 所以以下的算法只是用于演示数值解常微分方程的大致原理. 这种方法可以轻易地拓展到多个天体的情况, 而多体情况没有一般的解析解. 一种效率更高的常见常微分方程数值算法可参考 “四阶龙格库塔法”.

   直角坐标系中, 设中心天体质量为 $M$, 固定在原点不动.根据牛顿万有引力定律,质量为 $m$ 的行星受到中心天体的力为

\begin{equation} \boldsymbol{\mathbf{F}} = -G \frac{Mm}{r^2} \hat{\boldsymbol{\mathbf{r}}} = -G\frac{Mm}{r^3} \boldsymbol{\mathbf{r}} \end{equation}
其中 $ \boldsymbol{\mathbf{r}} $ 为行星的位矢(设行星在 $xy$ 平面上运动). 根据牛顿第二定律, 加速度为
\begin{equation} \boldsymbol{\mathbf{a}} = \frac{ \boldsymbol{\mathbf{F}} }{m} = -G\frac{M}{r^3} \boldsymbol{\mathbf{r}} \end{equation}
用 $r = \sqrt{x^2+y^2}$ 以及 $ \boldsymbol{\mathbf{r}} = x \hat{\boldsymbol{\mathbf{x}}} + y \hat{\boldsymbol{\mathbf{y}}} $, 其中 $x,y$ 看成 $t$ 的函数. 考虑到 $a_x = \ddot x$, $a_y = \ddot y$, 可以列出二阶微分方程组(变量上方两点表示关于时间的二阶导数)
\begin{equation} \left\{\begin{aligned} \ddot{x} = -\frac{GM}{(x^2+y^2)^{3/2}}x\\ \ddot{y} = -\frac{GM}{(x^2+y^2)^{3/2}}y \end{aligned}\right. \end{equation}
假设已知初值条件 $x(0) = x_0$, $y(0) = y_0$, $\dot x(0) = v_{x0}$, $\dot y(0) = v_{y0}$. 下面用 “弹簧振子受迫振动的简单数值计算” 中类似的方法求接下来行星的运动轨迹.

   1.将初始条件代入式 3 ,得到初始加速度

\begin{equation} \left\{\begin{aligned} \ddot x(0) &= -\frac{GM}{(x_0^2 + y_0^2)^{3/2}} x_0\\ \ddot y(0) &= -\frac{GM}{(x_0^2 + y_0^2)^{3/2}} y_0 \end{aligned}\right. \end{equation}

   2.设经过一段极微小的时间步长 $\Delta t$ (例如 $0.0001$, 数值越小误差越小), 根据微分近似(微分近似在这里的物理意义是在 $\Delta t$ 内速度和加速度都近似为常数)

\begin{equation} \begin{cases} \dot x(\Delta t) \approx \ddot x(0)\Delta t + x(0)\\ \dot y(\Delta t) \approx \ddot y(0)\Delta t + \dot y(0) \end{cases} \qquad \begin{cases} x(\Delta t) \approx \dot x(0)\Delta t + x(0)\\ y(\Delta t) \approx \dot y(0)\Delta t + y(0) \end{cases} \end{equation}

   3.把 $x(\Delta t), y(\Delta t)$ 再次代入式 3 , 得到 $x''(\Delta t), y''(\Delta t)$, 再次利用微分近似求出 $x(2\Delta t), y(2\Delta t), x'(2\Delta t), y'(2\Delta t) \dots$ 如此循环下去就可以得到每隔 $\Delta t$ 的数值解.

代码 1:kepler.m
% 参数设定
GM = 1; % 万有引力常数乘以中心天体质量
x0 = 1; y0 = 0; % 初始位置
vx0 = 0; vy0 = 0.7; % 初始速度
T = 4; Nstep = 4000; % 总时间和步数
dt = T/Nstep; % 步长

% 矩阵预赋值
x = nan(Nstep,1); y = x;
x1 = x;  y1 = x;
x2 = x; y2 = x;

% 初始位置,速度,加速度
x(1) = x0; y(1) = y0; % 初位置
x1(1) = vx0; y1(1) = vy0; % 初速度
x2(1) = -GM*x(1)/(x(1)^2+y(1)^2)^(3/2); % 代入方程得到 x''(0)
y2(1) = -GM*y(1)/(x(1)^2+y(1)^2)^(3/2); % 代入方程得到 y''(0)

% 迭代循环
for ii = 2:Nstep
    x(ii) = x(ii-1)+x1(ii-1)*dt; % x的微分
    y(ii) = y(ii-1)+y1(ii-1)*dt; % y的微分
    
    x1(ii) = x1(ii-1)+x2(ii-1)*dt; % x' 的微分
    y1(ii) = y1(ii-1)+y2(ii-1)*dt; % y' 的微分

    x2(ii) = -GM*x(ii)/(x(ii)^2+y(ii)^2)^(3/2); % 代入微分方程求出 x''
    y2(ii) = -GM*y(ii)/(x(ii)^2+y(ii)^2)^(3/2); % 代入微分方程求出 y''
end

% 画图
plot(x,y); % 画行星轨道
axis equal; % xy坐标长度一致
hold on; % 继续画图
scatter(0,0); % 标出中心天体

   程序运行结果如图 1 所示. 注意行星轨道并不是一个闭合的椭圆, 这是由于这种算法误差较大, 为了减小误差, 可以增加程序中 Nstep 的值.

图
图 1:运行结果
致读者: 小时物理百科一直以来坚持所有内容免费无广告,这导致我们处于日渐严重的亏损状态。长此以往很可能会最终导致我们不得不选择商业化,例如大量广告,内容付费, 甚至被收购。 因此,我们鼓起勇气在此请求广大读者热心捐款,使网站得以健康发展。如果看到这条信息的每位读者能慷慨捐助 10 元,我们一个星期内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。 但遗憾的是只有不到 1% 的读者愿意捐款, 他们的付出帮助了 99% 的读者免费获得知识, 我们在此表示感谢。

         

© 小时物理百科 保留一切权利