MATLAB教学：模拟三体运动

代码如下

clc;clear;
figure(‘name’,’三体运动’);%设置标题名字
N=3;
e=2.63;%电荷
%x=zeros(1,N);y=zeros(1,N);vx=zeros(1,N);vy=zeros(1,N);ax=zeros(1,N);ay=zeros(1,N);ke=ones(1,N);
x=[-sqrt(3),0,sqrt(3)];
y=[-1,2,-1];%设置三个质点的初始位置
vx=[-0.5,1,-0.5];
vy=[sqrt(3)/2,0,-sqrt(3)/2];%设置三个质点的初始速度
ax=zeros(1,N);
ay=zeros(1,N);
ke=[e,e,e];%设置三个质点的电荷相对值
M=[2,2,2];%设置三个质点的质量相对值
dt=0.005;
pausetime=0.002;%设置时间微小元长度，越小演示越精细，但越慢；设置暂停时间；

set(gcf,’doublebuffer’,’on’)%消除抖动
set(gca,’xlim’,[-3 3],’ylim’,[-3 3]);%设置坐标轴范围
hold on;
axis equal;
for m=1:N
p(m)=plot(x(m),y(m),’color’,’b’,’marker’,’.’,’markersize’,15); %所有质点初始位置以及大小设置
end
for jj=1:50000 %设定运行距离
for m=1:N
ax(m)=0;ay(m)=0;
for n=1:N
if m~=n
ax(m)=ax(m)+ke(n)*ke(m)*(x(n)-x(m))*((x(n)-x(m))^2+(y(n)-y(m))^2)^(-1.5)/M(m);
%按吸引力的格式写的加速度，如果要改为排斥力，需要将等号后面的m和n交换位置
ay(m)=ay(m)+ke(n)*ke(m)*(y(n)-y(m))*((x(n)-x(m))^2+(y(n)-y(m))^2)^(-1.5)/M(m);
else
end
end
x(m)=x(m)+vx(m)*dt+0.5*ax(m)*dt^2;%计算质点的新位置
y(m)=y(m)+vy(m)*dt+0.5*ay(m)*dt^2;
vx(m)=vx(m)+ax(m)*dt;
vy(m)=vy(m)+ay(m)*dt;
set(p(m),’xdata’,x(m),’ydata’,y(m));%设置质点的运动过程
if m==1
plot(x(m),y(m),’b*’);%画出三个质点的运动轨迹
elseif m==2
plot(x(m),y(m),’r*’);
else
plot(x(m),y(m),’g*’);
end
if abs(x(m))>1000||abs(y(m))>1000%如果质点已经运动到边框外面则停止运行，跳出该层循环
break;
end
end
if abs(x(m))>1000||abs(y(m))>1000 %如果质点已经运动到边框外面则停止运行，停止运行
break;
end
pause(pausetime);%暂停一会
drawnow
end