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

蒲丰氏投针问题的模拟过程，随机数发生器也是自编的，以供大家参考和提出建议。谢谢。（seed1和seed2最好选择3和5，为了使投针次数达到1000000，CVF进行如下设置Project-settings-link- output，将stack allocations reserve:设为1000000000） program getpi implicit none real,parameter:a=5,L=4,pi=3.14159 integer:n1,i,counter=0 real,allocatable:R1(:),R2(:) real:theta,x,pi1 write(*,*) input the size of the array: read(*,*) n1 allocate(R1(n1) allocate(R2(n1) call random(n1,R1,R2) do i=1,n1 x=a*(2*R1(i)-1) theta=pi*R2(i) if(abs(x)L*sin(theta) counter=counter+1 end do pi1=(1.0*n1)/(counter*1.0)*(2.0*L/a) write(*,*) this is PI: write(*,*) pi write(*,(this is ratio of counter to total number,F10.6) 1.0 &*counter/n1 stop end program subroutine random(n,R1,R2) implicit none integer n,seed1,seed2,i,little,m real:R1(n),R2(n) integer:temp1(n),temp2(n) write(*,*) input seed1 read(*,*) seed1 write(*,*) input seed2 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(R1(1)0.5) little=little+1 do i=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 . end do ; write(*,*) ratio of number which is little than 0.5 write(*,*) 1.0*little/n return end subroutine拟蒙特卡洛方法(Quasi-Monte Carlo)积分实例%使用Matlab提供的函数求积分,exp(-1/2*x2)在(0,1)间积分 format long;syms xa = sym(1/2);f = exp(-a*x2);ezplot(f)disp(int(f,-1,1);fprintf(integral result:%1.18f.n,double(int(f,0,1);%disp(double(int(f,0,1);复制代码%使用拟蒙特卡洛方法积分%得到拟蒙特卡洛序列,即低偏差序列,halton法%如果有相关的工具箱的话,可以用Matlab里面的haltonset,faureset,sobolset函数实现, x=halton(10000,2,5577);n=length(x);mju=0;for i=1:n mju=mju + exp(-0.5*x(i)2);endmju=mju/n;fprintf(Quasi-Monte Carlo result:%1.18f.n,mju);%disp(mju);%使用蒙特卡洛方法积分%得到Uniform序列,x=random(unif,0,1,10000,1);n=length(x);mju=0;for i=1:n mju=mju + exp(-0.5*x(i)2);endmju=mju/n;fprintf(Monte Carlo result:%1.18f.n,mju);%=生成HALTON序列=function result = halton( m,base,seeder )%生成HALTON序列% Check inputsif nargin 3seeder = 0;if nargin 0 element(k)=rem(temp,base(j); temp=fix(temp/base(j); k=k