三體運(yùn)動(dòng)的matlab演示

1、-作者xxxx-日期xxxx三體運(yùn)動(dòng)的matlab演示【精品文檔】三體運(yùn)動(dòng)的matlab演示figure(name,三體運(yùn)動(dòng));%設(shè)置標(biāo)題名字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;%設(shè)置三個(gè)質(zhì)點(diǎn)的初始位置vx=-1.5,1,-0.8;vy=1.2,1,0.6;%設(shè)置三個(gè)質(zhì)點(diǎn)的初始速度ax=zeros(1,N);ay=zeros(1,N);ke=1.5,6,2;%設(shè)置三個(gè)質(zhì)點(diǎn)的電荷相對(duì)值M=1,5,1;

2、%設(shè)置三個(gè)質(zhì)點(diǎn)的質(zhì)量相對(duì)值dt=0.005;pausetime=0.002;%設(shè)置時(shí)間微小元長(zhǎng)度，越小演示越精細(xì)，但越慢；設(shè)置暫停時(shí)間； set(gcf,doublebuffer,on) %消除抖動(dòng)set(gca,xlim,-7 7,ylim,-7 7);%設(shè)置坐標(biāo)軸范圍hold on;axis equal;for m=1:N p(m)=plot(x(m),y(m),color,k,marker,.,markersize,15); %所有質(zhì)點(diǎn)初始位置以及大小設(shè)置endfor jj=1:5000 %設(shè)定運(yùn)行距離for m=1:N ax(m)=0;ay(m)=0; for n=1:N if m=n

3、 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);%按吸引力的格式寫(xiě)的加速度，如果要改為排斥力，需要將等號(hào)后面的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;%計(jì)算質(zhì)點(diǎn)的新位置 y(m)=y(m)+vy(m)*dt+0.5*ay(m)*dt2; vx(m)=vx(m)+ax(m)*dt; vy(m)

4、=vy(m)+ay(m)*dt; set(p(m),xdata,x(m),ydata,y(m);%設(shè)置質(zhì)點(diǎn)的運(yùn)動(dòng)過(guò)程 plot(x(m),y(m),color,b);%畫(huà)出三個(gè)質(zhì)點(diǎn)的運(yùn)動(dòng)軌跡 if abs(x(m)10|abs(y(m)10 %如果質(zhì)點(diǎn)已經(jīng)運(yùn)動(dòng)到邊框外面則停止運(yùn)行，跳出該層循環(huán) break; endendif abs(x(m)10|abs(y(m)10 %如果質(zhì)點(diǎn)已經(jīng)運(yùn)動(dòng)到邊框外面則停止運(yùn)行，停止運(yùn)行 break;end % pause(pausetime); %暫停一會(huì)drawnowend運(yùn)行結(jié)果說(shuō)明：更改不同的參數(shù)得到不同的運(yùn)行結(jié)果。以上程序只是平面內(nèi)的演示，希望讀者根據(jù)平面內(nèi)的模擬程序得到空間內(nèi)的演示程序。也可以添加更多的質(zhì)點(diǎn)，得到更多體的運(yùn)動(dòng)，不過(guò)會(huì)使得運(yùn)行變慢，較好的計(jì)算機(jī)才能做到。如要添加到四個(gè)質(zhì)點(diǎn)，則需要N=4;且x,y,vx,vy,ke,M都要相應(yīng)的有四個(gè)初始值。程序目前存在的問(wèn)題：當(dāng)兩個(gè)質(zhì)點(diǎn)運(yùn)動(dòng)到一點(diǎn)（即發(fā)生碰