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

clear

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

close all

% 三個(gè)恒星的質(zhì)量都是1

ms = 1 ;

mt = 1 ;

mj = 1 ;

% 無(wú)量綱后萬(wàn)有引力常數(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ù)是加速，牛頓第二定律

% 以及萬(wàn)有引力定律得到常微分方程組

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求解常微分方程組的初值問(wèn)題

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

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

???????????????Y(i,1),Y(i,2),'-or','LineWidth',4, ... %畫(huà)第一個(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")