三體系統(tǒng)運(yùn)動(dòng)軌跡腳本如下，感興趣的朋友不妨下載北太天元運(yùn)行試一試！

%模擬三個(gè)恒星組成的系統(tǒng)的三體運(yùn)動(dòng)clear

load_plugin("time");  %為了使用北太天元軟件的pause插件函數(shù)

close all

% 三個(gè)恒星的質(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];

%初始時(shí)刻

to = 0;

%計(jì)算終止時(shí)刻

tf = 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]))); %三個(gè)質(zhì)點(diǎn)的x坐標(biāo)（在所有時(shí)刻)的最小值

xmax = max(max(Y(:,[1,3,5])));

ymin = min(min(Y(:,[2,4,6]))); %三個(gè)質(zhì)點(diǎn)的y坐標(biāo)（在所有時(shí)刻)的最小值

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', ... %畫一個(gè)白色的斜線代替axis([xmin,xmax,ymin,ymax])設(shè)置畫圖范圍

               Y(i-shift:i,1),Y(i-shift:i,2),'r','LineWidth',2, ... %畫第一個(gè)恒星在i-shift個(gè)時(shí)刻和第i個(gè)時(shí)刻件的軌跡

               Y(i,1),Y(i,2),'-or','LineWidth',4, ... %畫第一個(gè)恒星在第i個(gè)時(shí)刻所在的位置

               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('時(shí)間=%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")