matlab模拟三体运动_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

打赏作者

H.C

微信支付

支付宝