三体运动的matlab演示

figure('name','三体运动');%设置标题名字

N=3; %x=zeros(1,N);y=zeros(1,N);vx=zeros(1,N);vy=zeros(1,N);ax=zer os(1,N);ay=zeros(1,N);ke=ones(1,N);

x=[0,2,3];y=[0,2,0];%设置三个质点的初始位置

vx=[-1.5,1,-0.8];vy=[1.2,1,0.6];%设置三个质点的初始速度

ax=zeros(1,N);ay=zeros(1,N);

ke=[1.5,6,2];%设置三个质点的电荷相对值

M=[1,5,1];%设置三个质点的质量相对值

dt=0.005;pausetime=0.002;%设置时间微小元长度，越小演示越精细，但越慢；设置暂停时间；

set(gcf,'doublebuffer','on') %消除抖动

set(gca,'xlim',[-7 7],'ylim',[-7 7]);%设置坐标轴范围

hold on;axis equal;

for m=1:N

p(m)=plot(x(m),y(m),'color','k','marker','.','markersize',15); %所有质点初始位置以及大小设置

end

for jj=1:5000 %设定运行距离

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

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));%设置质点的运动过程

plot(x(m),y(m),'color','b');%画出三个质点的运动轨迹

if abs(x(m))>10||abs(y(m))>10 %如果质点已经运动到边框外面则停止运行，跳出该层循环

break;

end

if abs(x(m))>10||abs(y(m))>10 %如果质点已经运动到边框外面则停止运行，停止运行break;

end

% pause(pausetime); %暂停一会

drawnow

end

运行结果

说明：更改不同的参数得到不同的运行结果。

以上程序只是平面内的演示，希望读者根据平面内的模拟程序得到空间内的演示程序。

也可以添加更多的质点，得到更多体的运动，不过会使得运行变慢，较好的计算机才能做到。如要添加到四个质点，则需要N=4;且x,y,vx,vy,ke,M都要相应的有四个初始值。

程序目前存在的问题：

当两个质点运动到一点（即发生碰撞时），会产生速度的突变，如何能较好解决这个问题？如果能够，希望能得到空间内的演示。

见文库内本人所编写其他的matlab模拟。