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

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

2、%設(shè)置三個質(zhì)點(diǎn)的質(zhì)量相對值dt=0.005;pausetime=0.002;%設(shè)置時間微小元長度，越小演示越精細(xì)，但越慢；設(shè)置暫停時間； set(gcf,doublebuffer,on) %消除抖動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);%按吸引力的格式寫的加速度，如果要改為排斥力，需要將等號后面的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;%計算質(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)動過程 plot(x(m),y(m),color,b);%畫出三個質(zhì)點(diǎn)的運(yùn)動軌跡 if abs(x(m)10|abs(y(m)10 %如果質(zhì)點(diǎn)已經(jīng)運(yùn)動到邊框外面則停止運(yùn)行，跳出該層循環(huán) break; endendif abs(x(m)10|abs(y(m)10 %如果質(zhì)點(diǎn)已經(jīng)運(yùn)動到邊框外面則停止運(yùn)行，停止運(yùn)行 break;end % pause(pausetime); %暫停一會drawnowend運(yùn)行結(jié)果說明：更改不同的參數(shù)得到不同的運(yùn)行結(jié)果。以上程序只是平面內(nèi)的演示，希望讀者根據(jù)平面內(nèi)的模擬程序得到空間內(nèi)的演示程序。也可以添加更多的質(zhì)點(diǎn)，得到更多體的運(yùn)動，不過會使得運(yùn)行變慢，較好的計算機(jī)才能做到。如要添加到四個質(zhì)點(diǎn)，則需要N=4;且x,y,vx,vy,ke,M都要相應(yīng)的有四個初始值。程序目前存在的問題：當(dāng)兩個質(zhì)點(diǎn)運(yùn)動到一點(diǎn)（即發(fā)生碰