蒙特卡罗方法C.ppt

clearall;N=1000000;n=0;
fori=1:N
x1=rand;x2=rand;x3=rand;x4=rand;
r2=x1^2+x2^2+x3^2+x4^2;
ifr2<=1
n=n+1;
end
end
P=n/N;SigmaP=sqrt(P*(1-P))/sqrt(N);
V=P*2^4
SigmaV=SigmaP*V
V/pi/pi
SigmaV/pi/pi§5.2一维定积分计算的平均值法
§5.3半解析方法
§5.4重要抽样法
§5.5控制变量法
§5.6拟蒙特卡罗方法
超均匀随机数
范德科普特（VanderCorput）随机数产生器clearall;N=10000;n=0;Sumx=0;Sumxx=0;
fori=1:N
x1=rand;x2=rand;
ifx2<=x1
n=n+1;
end
Sumx=Sumx+x1;
Sumxx=Sumxx+x1*x1;
end
P=n/N;Ipoint=P
SigmaP=sqrt(P*(1-P))/sqrt(N);SigmaIpoint=SigmaP
Iav=Sumx/N
SigmaIav=sqrt(Sumxx/N-Iav*Iav)/sqrt(N)十进制整数§5.7分层抽样
简单分层抽样
最优分层抽样clearall;N=10000;k=10;m=N/k;deltx=1.0/k;
Sumx(1:k)=0;Sumxx(1:k)=0;
fori=1:k
forj=1:m
x=(i-1)*deltx+rand*deltx;
Sumx(i)=Sumx(i)+x;
Sumxx(i)=Sumxx(i)+x*x;
end
end
Ik=Sumx/m
SigmaIk=sqrt(Sumxx/m-Ik.*Ik)/sqrt(m)
Iav=sum(Ik)/k
SigmaIav=sqrt(sum(SigmaIk.*SigmaIk))/kclearall;N=10000;k=5;m=N/k;deltx=1.0/k;
Sumf(1:k)=0;Sumff(1:k)=0;
fori=1:k
forj=1:m
x=(i-1)*deltx+rand*deltx;f=x^4;
Sumf(i)=Sumf(i)+f;
Sumff(i)=Sumff(i)+f*f;
end
end
Ik=Sumf/m
SigmaIk=sqrt(Sumff/m-Ik.*Ik)/sqrt(m)
Iav=sum(Ik)/k
SigmaIav=sqrt(sum(SigmaIk.*SigmaIk))/k第六章分子物理学中的模拟实验设两维容器中有N个粒子，初始速率相等。
粒子与容器壁碰撞时，速度只改变方向，不改变大小（即粒子与容器壁无能量交换）。
粒子与粒子碰撞时，为弹性碰撞。（1） 设定N个粒子的初始速率；
（2） 随机选择两个粒子（粒子j和粒子k），两粒子运动方向随机；
（3） 从实验室系转换到质心系，描述入射粒子：（4）质心系中两粒子进行碰撞：clearall;
N=10000;v0=1;
v(1:N)=v0;
deltv=3*v0/20;
fori=1:20
vbin(i)=(i-0.5)*deltv;
end
subplot(3,1,1),hist(v,vbin);
xlabel('v');ylabel('N');
fori=1:10*N
j=ceil(rand*N);k=ceil(rand*N);
whilek==j
k=ceil(rand*N);
end
phij=2*pi*rand;phik=2*pi*rand;
vj=[v(j)*cos(phij);v(j)*sin(phij)];
vk=[v(k)*cos(phik);v(k)*sin(phik)];
vc=1/2*(vj+vk);
wj=vj-vc;wk=vk-vc;10000个粒子在（1）初始时刻；（2）碰撞10000次后；（3）碰撞100000次后的速率分布质量不同的两种粒子的碰撞clearall;
N=100000;v0=1;m1=1;m2=25;r=4/5;
v(1:N)=v0;m(1:N*r)=m1;m(N*r+1:N)=m2;
deltv=0.1;deltE=0.1;
fori=1:80
vbin(i)=(i-0.5)*deltv;
Ebin(i)=(i-0.5)*deltE;
end
subplot(2,2,1),hist(v,vbin);
xlabel('v');ylabel('N');
subplot(2,2,2),hist((1/2*m.*v.^2).^(1/2),Ebin);
xlabel('sqrt(E)');ylabel('N');
fori=1:50*N
j=ceil(rand*N);k=ceil(rand*N);
whilek==j