图

天体运动的简单数值计算

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

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

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

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

\begin{equation} \leftgroup{ \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{equation}

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

\begin{equation} \leftgroup{ \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)} \qquad \leftgroup{ x(\Delta t) &\approx \dot x(0)\Delta t + x(0)\\ y(\Delta t) &\approx \dot y(0)\Delta t + y(0)} \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$ 的数值解.

  kepler.m

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
% 参数设定
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 元,我们几天内就能脱离亏损状态,并保证网站能在接下来的一整年里向所有读者继续免费提供优质内容。感谢您的支持。
—— 小时(项目创始人)

编辑词条 返回目录 返回主页 捐助项目 © 小时物理百科 保留一切权利