三体系统运动轨迹

三体系统运动轨迹脚本如下，感兴趣的朋友不妨下载北太天元运行试一试！

%模拟三个恒星组成的系统的三体运动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;

%计算终止时刻

tf = 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")