闲的蛋疼，学学matlab动画，三体运动还是有点问题（提示line没xdate属性，脑瓜也疼了），不多废话上代码

clc;clear;close;
figure('name','三体运动');
axis equal
grid on
hold on
colordef black
view(3)
xlabel('x')
ylabel('y')
zlabel('z')
hold on
%初始化参数设置(质量、初始位置、初始速度、万有引力常数、时间步长)
m1=rand()*10^23;
m2=rand()*10^23;
m3=rand()*10^23;

pos1=[rand()*10^7,rand()*10^7,rand()*10^7];
pos2=[rand()*10^7,rand()*10^7,rand()*10^7];
pos3=[rand()*10^7,rand()*10^7,rand()*10^7];
 
v1=[rand()*1000,rand()*1000,rand()*1000];
v2=[rand()*1000,rand()*1000,rand()*1000];
v3=[rand()*1000,rand()*1000,rand()*1000];
G=6.67*10^(-11);
dt=0.005;

%行星和轨迹（定义行星位置和轨迹变量）
planet1=plot3(pos1(1),pos1(2),pos1(3),'b:.','markersize',20);
planet2=plot3(pos2(1),pos2(2),pos2(3),'r:.','markersize',20);
planet3=plot3(pos3(1),pos3(2),pos3(3),'y:.','markersize',20);

h1=animatedline('color','b');
h2=animatedline('color','r');
h3=animatedline('color','y');
frame=0;
while true
     %不断循环计算距离、相互作用力、速度、加速度、位置
    r12=normest(pos2-pos1);
    r23=normest(pos3-pos2);
    r13=normest(pos3-pos1);
    
    F12_val=G*m1*m2/(r12^2);
    F23_val=G*m2*m3/(r23^2);
    F13_val=G*m1*m3/(r13^2);
    
    F12_dir=(pos2-pos1)/normest(pos2-pos1);
    F23_dir=(pos3-pos2)/normest(pos3-pos2);
    F13_dir=(pos3-pos1)/normest(pos3-pos1);
    
    a1=((F12_val)*(F12_dir)+(F13_val)*(F13_dir))/m1;
    a2=((-F12_val)*(F12_dir)+(-F23_val)*(F23_dir))/m2;
    a3=((-F13_val)*(F13_dir)+(-F23_val)*(F23_dir))/m3;
    
    v1=v1+a1*dt;
    v2=v2+a2*dt;
    v3=v3+a3*dt;
    
    pos1=pos1+v1*dt+1/2*a1*dt^2;
    pos2=pos2+v2*dt+1/2*a2*dt^2;
    pos3=pos3+v3*dt+1/2*a3*dt^2;
    
    frame=frame+1;
    if frame==10000
        frame=0;
        
       % set(planet1,'position',[pos1(1),pos1(2),pos1(3)]);
       
        %set(planet3,'position',pos3(1),'YData',pos3(2),'ZData',pos3(3));
        
        addpoints(h1,pos1(1),pos1(2),pos1(3));
        addpoints(h2,pos2(1),pos2(2),pos2(3));
        addpoints(h3,pos3(1),pos3(2),pos3(3));
        drawnow
    end
    if r12==0||r23==0||r13==0       %判断是否相撞
        break
    end
end

好吧很丑