（文章靈感來自盧朓老師的B站視頻）

最近電視劇《三體》的大熱，引起了大家對三體系統(tǒng)的注意力，今天就讓我們在北太天元上面模擬一下三體系統(tǒng)的運(yùn)動軌跡

首先，什么是三體系統(tǒng)呢？

三體（three-body problem）
天體力學(xué)中的基本力學(xué)模型。研究三個可視為質(zhì)點(diǎn)的天體在相互之間萬有引力作用下的運(yùn)動規(guī)律問題。
這三個天體的質(zhì)量、初始位置和初始速度都是任意的。
----------------------------摘自百度百科

簡單的說，就是我們需要模擬三個恒星組成的系統(tǒng)的三體運(yùn)動。

詳細(xì)代碼如下，代碼摘自互聯(lián)網(wǎng)

%模擬三個恒星組成的系統(tǒng)的三體運(yùn)動
clear
load_plugin("time"); 
%為了使用北太天元軟件的pause插件函數(shù)
close all
% 三個恒星的質(zhì)量都是1
ms = 1 ;
mt = 1 ;
mj = 1 ;
% 無量綱后萬有引力常數(shù)設(shè)置為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; 
%由位置的導(dǎo)數(shù)速度，速度的導(dǎo)數(shù)是加速，牛頓第二定律
% 以及萬有引力定律得到常微分方程組
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]))); 
%三個質(zhì)點(diǎn)的x坐標(biāo)（在所有時刻)的最小值
xmax = max(max(Y(:,[1,3,5])));
ymin = min(min(Y(:,[2,4,6]))); 
%三個質(zhì)點(diǎn)的y坐標(biāo)（在所有時刻)的最小值
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])設(shè)置畫圖范圍
            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")