三体运动的matlab演示.doc

三体运动的matlab演示figure(name,三体运动);%设置标题名字N=3; %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=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); %所有质点初始位置以及大小设置endfor 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 end x(m)=x(m)+vx(m)*dt+0.5*ax(m)*dt2;%计算质点的新位置 y(m)=y(m)+vy(m)*dt+0.5*ay(m)*dt2; 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; endendif abs(x(m)10|abs(y(m)10 %如果质点已经运动到边框外面则停止运行，停止运行 break;end % pause(pausetime); %暂停一会drawnowend运行结果说明：更改不同的参数得到不同的运行结果。以上程序只是平面内的演示，希望读者根据平面内的模拟程序得到空间内的演示程序。也可以添加更多的质点，得到更多体的运动，不过会使得运行变慢，较好的计算机才能做到。如要添加到四个质点，则需要N=4;且x,y,vx,vy,ke,M都要相应的有四个初始值。程序目前存在的问题：当两个质点运动到一点（