matlab&simulink仿真天体运行轨迹

1. 介绍

本博客来自北理工《计算机仿真与matlab》课程的大作业，有任何问题欢迎批评指正。
源代码请参见如下链接：
链接：https://pan.baidu.com/s/1hkAoi9c16_IoUYJSxX92hg
提取码：4gr6

2. 问题背景

我使用matlab和simulink建立天体系统的模型，并使用不同的数值积分方法仿真求解，模拟推演天体运动轨迹。

下文第3、4节将演示双星系统的仿真模拟，第5节将演示三体系统的仿真模拟。

3. 数学模型

行星所受的万有引力为：

$$\vec{F}=\frac{GMm}{{\left|\vec{R}-\vec{r}\right|}^2}\cdot \frac{\vec{R}-\vec{r}}{\left|\vec{R}-\vec{r}\right|}$$

其中 $G$、$M$、$m$ 分别为万有引力常量、恒星质量、行星质量，$\vec{R}$ 是恒星位置，$\vec{r}$ 是行星位置。

由牛二律建立行星运动微分方程：

$$\frac{d^2\vec{r}}{d{t}^2}=\frac{GM(\vec{R}-\vec{r})}{{\left|\vec{R}-\vec{r}\right|}^3}$$

考虑到恒星的质量远大于行星的质量，假设恒星在系统中静止，即$\vec{R}$ 为常矢量。因此上式为关于行星位置$\vec{r}$ 的二阶微分方程，求解时将$\vec{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
$\vec{R}\left(t\right)$= [0,0,0]
$\vec{r}\left(t\right)$ = [8e6,0,0]
$\vec{V}\left(t\right)$ = [0,0,0]
$\vec{v}\left(t\right)$ = [-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 数学模型

三体系统的数学模型中除新增一颗恒星外，与前文双星系统完全相同。以行星为例，行星的加速度为：
$$\vec{a}=\frac{GM\left(\overrightarrow{R_1}-\vec{r}\right)}{{\left|\overrightarrow{R_1}-\vec{r}\right|}^3}+\frac{GM\left(\overrightarrow{R_2}-\vec{r}\right)}{{\left|\overrightarrow{R_2}-\vec{r}\right|}^3}$$
同样可以建立系统的微分方程和状态空间方程，此处略。

5.2 仿真模拟

simulink模型图如下：

运行simulink，使用龙格-库塔法求解。参数及初始条件如下：

G = 6.67e-11
h = 0.001
M_1 = 25e6
M_2 = 50e6
m = 0.5e6
$\overrightarrow{R_1}\left(t\right)$ = [0,0,0]
$\overrightarrow{R_2}\left(t\right)$ = [0,100,200]
$\vec{r}\left(t\right)$ = [0,-100,300]
$\overrightarrow{V_1}\left(t\right)$ = [0,0,0]
$\overrightarrow{V_2}\left(t\right)$ = [400,0,0]
$\vec{v}\left(t\right)$ = [20,100,-50]

仿真模拟得到的运动轨迹如下：

图中蓝色轨迹为行星，红色轨迹为恒星1，黄色轨迹为恒星2。行星在漫漫长夜中逐渐被两颗恒星吸引，在围绕恒星2做短暂的圆周运动之后被甩出系统，重新进入漫漫长夜。