资源描述：
《蒲丰氏投针问题的模拟过程》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。
1、蒲丰氏投针问题的模拟过程，随机数发生器也是自编的，以供大家参考和提出建议。谢谢。（seed1和seed2最好选择3和5，为了使投针次数达到1000000，CVF进行如下设置Project->settings->link->output，将stackallocationsreserve:设为1000000000）programgetpi implicitnone real,parameter::a=5,L=4,pi=3.14159 integer::n1,i,counter=0
2、 real,allocatable::R1(:),R2(:) real::theta,x,pi1 write(*,*)'inputthesizeofthearray:' read(*,*)n1 allocate(R1(n1)) allocate(R2(n1)) callrandom(n1,R1,R2) doi=1,n1 x=a*(2*R1(i)-1) theta=pi*R2(i) if(abs(x)3、))counter=counter+1 enddo pi1=(1.0*n1)/(counter*1.0)*(2.0*L/a) write(*,*)'thisisPI:' write(*,*)pi write(*,"('thisisratioofcountertototalnumber',F10.6)")1.0 &*counter/n1 stop endprogram subroutinerandom(n,R1,R2) impli
4、citnone integern,seed1,seed2,i,little,m real::R1(n),R2(n) integer::temp1(n),temp2(n) write(*,*)'inputseed1' read(*,*)seed1 write(*,*)'inputseed2' read(*,*)seed2 m=2**30 m=m*2-1 temp1(1)=397204094 temp2(1)=
5、temp1(1) R1(1)=(1.0*temp1(1))/(1.0*m) R2(1)=(1.0*temp2(1))/(1.0*m) little=0 if(R1(1)<0.5)little=little+1 doi=1,n-1 temp1(i+1)=mod(seed1*temp1(i),m) R1(i+1)=(1.0*temp1(i+1))/(1.0*m) temp2(i+1)=mod(seed2*temp2(i),m) R2(
6、i+1)=(1.0*temp2(i+1))/(1.0*m) if(R1(i+1)<0.5)little=little+1. enddo; write(*,*)'ratioofnumberwhichislittlethan0.5' write(*,*)1.0*little/n return endsubroutine拟蒙特卡洛方法(Quasi-MonteCarlo)积分实例%使用Matlab提供的函数求积分,exp(-1/2*x^2)在(0,1)间积分forma
7、tlong;symsxa=sym(1/2);f=exp(-a*x^2);ezplot(f)disp(int(f,-1,1));fprintf('integralresult:%1.18f.',double(int(f,0,1)));%disp(double(int(f,0,1)));复制代码%使用拟蒙特卡洛方法积分%得到拟蒙特卡洛序列,即低偏差序列,halton法%如果有相关的工具箱的话,可以用Matlab里面的haltonset,faureset,sobolset函数实现,x=halton(100
8、00,2,5577);n=length(x);mju=0;fori=1:n mju=mju+exp(-0.5*x(i)^2);endmju=mju/n;fprintf('Quasi-MonteCarloresult:%1.18f.',mju);%disp(mju);%使用蒙特卡洛方法积分%得到Uniform序列,x=random('unif',0,1,10000,1);n=length(x);mju=0;fori=1:n mju=m
