1. 介绍
本博客来自北理工《计算机仿真与matlab》课程的大作业,有任何问题欢迎批评指正。
源代码请参见如下链接:
链接:https://pan.baidu.com/s/1hkAoi9c16_IoUYJSxX92hg
提取码:4gr6
2. 问题背景
我使用matlab和simulink建立天体系统的模型,并使用不同的数值积分方法仿真求解,模拟推演天体运动轨迹。
下文第3、4节将演示双星系统的仿真模拟,第5节将演示三体系统的仿真模拟。
3. 数学模型
行星所受的万有引力为:
F ⃗ = G M m ∣ R ⃗ − r ⃗ ∣ 2 ⋅ R ⃗ − r ⃗ ∣ R ⃗ − r ⃗ ∣ \vec{F}=\frac{GMm}{\left|\vec{R}-\vec{r}\right|^2}·\frac{\vec{R}-\vec{r}}{\left|\vec{R}-\vec{r}\right|} F=∣∣∣R−r∣∣∣2GMm⋅∣∣∣R−r∣∣∣R−r
其中 G G G、 M M M、 m m m 分别为万有引力常量、恒星质量、行星质量, R ⃗ \vec{R} R 是恒星位置, r ⃗ \vec{r} r 是行星位置。
由牛二律建立行星运动微分方程:
d 2 r ⃗ d t 2 = G M ( R ⃗ − r ⃗ ) ∣ R ⃗ − r ⃗ ∣ 3 \frac{d^2\vec{r}}{dt^2}=\frac{GM(\vec{R}-\vec{r})}{\left|\vec{R}-\vec{r}\right|^3} dt2d2r=∣∣∣R−r∣∣∣3GM(R−r)
考虑到恒星的质量远大于行星的质量,假设恒星在系统中静止,即 R ⃗ \vec{R} R 为常矢量。因此上式为关于行星位置 r ⃗ \vec{r} r 的二阶微分方程,求解时将 r ⃗ \vec{r} r 分解到x、y、z轴,得到与上式等价的三个微分方程组。
4. 仿真模拟
4.1 matlab仿真
函数pre_load_twobody.m,预定义仿真参数(仿真时间间隔h),天体参数(万有引力常量G和天体质量)和初始条件(天体的初始位置和初始速度)。
function [G,h,M,m,posM,posm, vM, vm] = pre_load_twobody()
% 宇宙参数
G = 6.67*10^-11; % 万有引力常量 N*m^2/kg^2
h = 0.001; % 时间元
% 天体的质量 kg
M = 1.9891*1e25;
m = 5.965*1e19;
% 天体的初始位置 矢量,m
posM = [0,0,0];
posm = [8*1e6,0,0];
% 天体的初始速度 矢量,m/s
vM = [0,0,0];
vm = [-2.5e3,7.5e3,0];
end
twobody.m,使用梯形法求解微分方程并绘制天体轨迹。(其中调用的count_a函数用于计算加速度,此处略,请参见源代码)
[G,h,M,m,posM,posm, vM, vm] = pre_load_twobody();
%% 设置绘制参数和绘制对象
colordef black
grid on
% 绘制天体
hold on
planetM = plot3(posM(1), posM(2), posM(3), 'r:.', 'markersize', 60);
planetm = plot3(posm(1), posm(2), posm(3), 'g:.', 'markersize', 40); % 梯形法计算的轨迹
% 绘制天体轨迹
hm = animatedline('color', 'g', 'LineWidth', 2);
legend([hm],{'梯形法'});
xlabel('X')
ylabel('Y')
zlabel('Z')
%% 作图循环
step = 0;
while true
[aM, am, r12] = count_a(G,M,m,posM,posm);
%预估下一时刻的位置。即由t时刻的x1、x2求x1p(t+1)
posM_p = posM + vM * h;
posm_p = posm + vm * h;
%预估下一个时刻的速度。即由t时刻的x1、x2求x2p(t+1)
vm_p = vm + am * h;
%计算am_p,用于梯形法更新速度
[aM_p, am_p, ~] = count_a(G,M,m,posM_p,posm_p);
% 更新下一时刻的位置。即由t时刻的x1、x2求x1(t+1)
posm = posm + (vm+vm_p)*h/2;
% 更新下一个时刻的速度。即由t时刻的x1、x2求x2(t+1)
vm = vm + (am+am_p)*h/2;
%%% 每个1000个frame绘图一下
step = step + 1;
if mod(step,1000)==0
% 通过set改变点的位置
set(planetM, 'Xdata', posM(1), 'Ydata', posM(2), 'Zdata', posM(3));
set(planetm, 'Xdata', posm(1), 'Ydata', posm(2), 'Zdata', posm(3));
addpoints(hm, posm(1), posm(2), posm(3));
view(37.5,30) % 查看三维的图。因为使用hold on后matlab会默认使用二维画布。
drawnow
end
% 若相撞就退出循环。为避免过零问题,判断相撞的条件应为距离小于某阈值。
if r12 <10^6
break
end
% 总仿真时间为2000s
if step*h==2000
break
end
end
运行twobody.m。运动轨迹随天体质量等参数设置和初始条件不同而不同。使用参数和初始条件如下:
G = 6.67e-11
h = 0.001
M = 2e25
m = 6e19
R
⃗
(
t
)
\vec{R}\left(t\right)
R(t)= [0,0,0]
r
⃗
(
t
)
\vec{r}\left(t\right)
r(t) = [8e6,0,0]
V
⃗
(
t
)
\vec{V}\left(t\right)
V(t) = [0,0,0]
v
⃗
(
t
)
\vec{v}\left(t\right)
v(t) = [-2.5e3,7.5e3,0];
轨迹绘制如下图,形成双星系统经典的椭圆轨道。
4.2. simulink仿真
simulink模型图如下:
simout模块用于输出行星轨迹以供后续画图。
gravitation模块用于计算行星的加速度,定义如下。
%% 此函数用于simulink,根据太阳位置和地球位置计算地球加速度
function a = gravitation(G,M,R,r)
% 计算距离
d = normest(R - r);
% 加速度
a_val = G * M / (d ^ 2);
a_dir = (R - r) / d;
a = a_val * (a_dir);
end
在Model Properties选项卡中(如下图)设置InitFcn回调函数,以在模型运行时自动将常量添加到工作区。
以从速度到位置的积分器模块为例,在该模块中设置初始条件为posm变量,该变量在InitFcn回调函数中定义。对于其他积分器初始条件的设置同理。
4.3 结果分析
运行simulink模型图,下图为取h = 0.001时的仿真结果。绿色轨迹为使用梯形法计算得到,黄色轨迹为使用龙格-库塔法得到。matlab的梯形法和simulink的龙格-库塔法所得到的行星轨迹几乎完全重叠。这说明尽管龙格-库塔法的精度高于梯形法,但当取的时间间隔很小时,梯形法也可以达到很好的效果,不过需要以较长的运行时间为代价。
为缩短仿真时间,下图为取h=0.1时的仿真结果。可以看到龙格-库塔法和梯形法出现了较明显的差别,说明时间间隔h的增加会导致截断误差的增加,导致梯形法的仿真效果明显不如龙格-库塔法好。
5. 三体系统
5.1 数学模型
三体系统的数学模型中除新增一颗恒星外,与前文双星系统完全相同。以行星为例,行星的加速度为:
a
⃗
=
G
M
(
R
1
⃗
−
r
⃗
)
∣
R
1
⃗
−
r
⃗
∣
3
+
G
M
(
R
2
⃗
−
r
⃗
)
∣
R
2
⃗
−
r
⃗
∣
3
\vec{a}=\frac{GM\left(\vec{R_1}-\vec{r}\right)}{\left|\vec{R_1}-\vec{r}\right|^3}\ +\ \frac{GM\left(\vec{R_2}-\vec{r}\right)}{\left|\vec{R_2}-\vec{r}\right|^3}\
a=∣∣∣R1−r∣∣∣3GM(R1−r) + ∣∣∣R2−r∣∣∣3GM(R2−r)
同样可以建立系统的微分方程和状态空间方程,此处略。
5.2 仿真模拟
simulink模型图如下:
运行simulink,使用龙格-库塔法求解。参数及初始条件如下:
G = 6.67e-11
h = 0.001
M_1 = 25e6
M_2 = 50e6
m = 0.5e6
R
1
⃗
(
t
)
\vec{R_1}\left(t\right)
R1(t) = [0,0,0]
R
2
⃗
(
t
)
\vec{R_2}\left(t\right)
R2(t) = [0,100,200]
r
⃗
(
t
)
\vec{r}\left(t\right)
r(t) = [0,-100,300]
V
1
⃗
(
t
)
\vec{V_1}\left(t\right)
V1(t) = [0,0,0]
V
2
⃗
(
t
)
\vec{V_2}\left(t\right)
V2(t) = [400,0,0]
v
⃗
(
t
)
\vec{v}\left(t\right)
v(t) = [20,100,-50]
仿真模拟得到的运动轨迹如下:
图中蓝色轨迹为行星,红色轨迹为恒星1,黄色轨迹为恒星2。行星在漫漫长夜中逐渐被两颗恒星吸引,在围绕恒星2做短暂的圆周运动之后被甩出系统,重新进入漫漫长夜。