三体系统在北太天元上的模拟

（文章灵感来自卢朓老师的B站视频）

最近电视剧《三体》的大热，引起了大家对三体系统的注意力，今天就让我们在北太天元上面模拟一下三体系统的运动轨迹

首先，什么是三体系统呢？

三体（three-body problem）
天体力学中的基本力学模型。研究三个可视为质点的天体在相互之间万有引力作用下的运动规律问题。
这三个天体的质量、初始位置和初始速度都是任意的。
----------------------------摘自百度百科

简单的说，就是我们需要模拟三个恒星组成的系统的三体运动。

详细代码如下，代码摘自互联网

%模拟三个恒星组成的系统的三体运动
clear
load_plugin("time"); 
%为了使用北太天元软件的pause插件函数
close all
% 三个恒星的质量都是1
ms = 1 ;
mt = 1 ;
mj = 1 ;
% 无量纲后万有引力常数设置为1
G = 1 ;
%初始条件 [xs,ys,xt,yt,xj,yj,vxs,vys,vxt,vyt,vjx,vjt]
CI = [0 -0.1 2 2 5 0 0 0 0 0 0 0];
%初始时刻to = 0;
%计算终止时刻t
f = 120; 
%由位置的导数速度，速度的导数是加速，牛顿第二定律
% 以及万有引力定律得到常微分方程组
fxy = @(ps, pt, pj,ms,mt,mj) ...G*( mt.*(pt-ps)./norm(pt-ps).^3 ...+ mj.*(pj-ps)./norm(pj-ps).^3 );
F = @(t,Y) [Y(7);Y(8);Y(9);Y(10);Y(11);Y(12);...
      fxy(Y([1,2]),Y([3,4]),Y([5,6]),ms,mt,mj); ...
      fxy(Y([3,4]),Y([1,2]),Y([5,6]),mt,ms,mj); ...
      fxy(Y([5,6]),Y([3,4]),Y([1,2]),mj,mt,ms); ...];
%使用ode45求解常微分方程组的初值问题
[t,Y]=ode45(F,[to,tf],CI); 
%plot(Y(:,1),Y(:,2),'r',Y(:,3),Y(:,4),'g',Y(:,5),Y(:,6),'b')
yo = Y(1) ;
dto = 0.3 ;
plotmax = 100 ;
T=to ; 
xmin = min(min(Y(:,[1,3,5]))); 
%三个质点的x坐标（在所有时刻)的最小值
xmax = max(max(Y(:,[1,3,5])));
ymin = min(min(Y(:,[2,4,6]))); 
%三个质点的y坐标（在所有时刻)的最小值
ymax = max(max(Y(:,[2,4,6]))); 
clf
close all
figure('Position',[0 0 1550 800])
hold off
told = 0;
for i = 1:length(Y(:,1))
    dt = abs(Y(i,1)-yo)/abs(Y(i,7));
    if dt >= dto
        if i>plotmax
            shift = plotmax;
        else
            shift = i-1;
        end
        plot(...
            [xmin,xmax],[ymin,ymax], 'w', ... %画一个白色的斜线代替axis([xmin,xmax,ymin,ymax])设置画图范围
            Y(i-shift:i,1),Y(i-shift:i,2),'r','LineWidth',2, ... %画第一个恒星在i-shift个时刻和第i个时刻件的轨迹
            Y(i,1),Y(i,2),'-or','LineWidth',4, ... %画第一个恒星在第i个时刻所在的位置
            Y(i-shift:i,3),Y(i-shift:i,4),'g','LineWidth',2, ...
            Y(i,3),Y(i,4),'-og','LineWidth',4, ...
            Y(i-shift:i,5),Y(i-shift:i,6),'b','LineWidth',2, ...
            Y(i,5),Y(i,6),'-ob','LineWidth',4)
        title(sprintf('时间=%f',t(i)))
        T=[T;t(i)];
        yo = Y(i,1) ;
        vo = Y(i,7) ;
        end
 pause(0.01)
 end
 X=[0:1:length(T)-1];
 figure(2)
 plot(X,T)
 plot(Y(:,1),Y(:,2),'r', 'LineWidth',2, ...
 Y(:,3),Y(:,4),'g','LineWidth',2, ...
 Y(:,5),Y(:,6),'b', 'LineWidth',2)
 unload_plugin("time")