蒲丰氏投针问题的模拟过程

《蒲丰氏投针问题的模拟过程》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

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      real,allocatable::

2、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、unter*1.0)*(2.0*L/a)    write(*,*)'thisisPI:'      write(*,*)pi    write(*,"('thisisratioofcountertototalnumber',F10.6)")1.0    &*counter/n1    stop      endprogram        subroutinerandom(n,R1,R2)    implicitnone      integern,seed1,seed2,i,little,m    real::R1(n),R2(n)      i

4、nteger::temp1(n),temp2(n)    write(*,*)'inputseed1'      read(*,*)seed1      write(*,*)'inputseed2'    read(*,*)seed2      m=2**30      m=m*2-1      temp1(1)=397204094      temp2(1)=temp1(1)      R1(1)=(1.0*temp1(1))/(1.0*m)      R2(1)=(1.0*temp2(1))/(1.0*m)    little=0    if(R

5、1(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(i+1)=(1.0*temp2(i+1))/(1.0*m)    if(R1(i+1)<0.5)little=little+1.      enddo;    write(*,*)'ratioofnumberwhichislittletha

6、n0.5'      write(*,*)1.0*little/n    return    endsubroutine拟蒙特卡洛方法(Quasi-MonteCarlo)积分实例%使用Matlab提供的函数求积分,exp(-1/2*x^2)在(0,1)间积分formatlong;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)));

7、复制代码%使用拟蒙特卡洛方法积分%得到拟蒙特卡洛序列,即低偏差序列,halton法%如果有相关的工具箱的话,可以用Matlab里面的haltonset,faureset,sobolset函数实现,x=halton(10000,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

8、('unif',0,1,10000,1);n=length(x);mju=0;fori=1:n  mju=m