基于MATLAB计算粒子速度分布

clear;clc N=10000;v0=1; %粒子个数，初始速度 v(1:N)=v0; %定义所有粒子的初始速度为v0 deltv=3*v0/20; for i=1:20 vbin(i)=(i-0.5)*deltv; %定义一个速度分布区间，用于绘图 end subplot(3,1,1),hist(v,vbin); xlabel('v') ylabel('N') title('初始速度分布') for i=1:10*N %碰撞总次数为10*N j=ceil(rand*N); %随机选取粒子j k=ceil(rand*N); %随机选取粒子k while k==j %如果粒子k和j相同，则重新选粒子k，直到k和j不是同一个粒子 k=ceil(rand*N); end phij=2*pi*rand; phik=2*pi*rand; vj=[v(j)*cos(phij);v(j)*sin(phij)]; %粒子j在x方向和y方向的初始速度 vk=[v(k)*cos(phik);v(k)*sin(phik)]; %粒子k在x方向和y方向的初始速度 vc=(vj+vk)/2; %粒子i和j组成的质心速度（粒子质量为1） wj=vj-vc; %粒子j在质心系里的初始速度 wk=vk-vc; %粒子k在质心系里的初始速度 phic=rand*2*pi; cosphic=cos(phic); sinphic=sin(phic); vj=[cosphic,-sinphic;sinphic,cosphic]*wj+vc; %实验室系下，粒子j碰撞之后的速度，包括vx分量和vy分量 vk=[cosphic,-sinphic;sinphic,cosphic]*wk+vc; %实验室系下，粒子k碰撞之后的速度，包括vx分量和vy分量 v(j)=sqrt(vj'*vj); %求出粒子j在实验室系下的速度（x方向速度和y方向速度合成） v(k)=sqrt(vk'*vk); %求出粒子k在实验室系下的速度（x方向速度和y方向速度合成） if i==N %输出碰撞N次之后的速度分布 subplot(3,1,2),hist(v,vbin); xlabel('v') ylabel('N') title([num2str(N) '次碰撞之后，速度分布']) end end subplot(3,1,3),hist(v,vbin); xlabel('v') ylabel('N') title([num2str(10*N) '次碰撞之后，速度分布']) %绘制与理论值比较的曲线 nv=hist(v,vbin); p=nv/N/deltv; for i=1:20 f(i)=exp(-vbin(i)^2/v0^2)*2*vbin(i)/v0^2; end figure plot(vbin,p,'o',vbin,f,'-') xlabel('v') ylabel('f(v)') title('理论曲线比较')