蒲豐氏投針問(wèn)題的模擬過(guò)程

《蒲豐氏投針問(wèn)題的模擬過(guò)程》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在行業(yè)資料-天天文庫(kù)。

1、蒲豐氏投針問(wèn)題的模擬過(guò)程，隨機(jī)數(shù)發(fā)生器也是自編的，以供大家參考和提出建議。謝謝。（seed1和seed2最好選擇3和5，為了使投針次數(shù)達(dá)到1000000，CVF進(jìn)行如下設(shè)置Project->settings->link->output，將stackallocationsreserve:設(shè)為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)積分實(shí)例%使用Matlab提供的函數(shù)求積分,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、復(fù)制代碼%使用擬蒙特卡洛方法積分%得到擬蒙特卡洛序列,即低偏差序列,halton法%如果有相關(guān)的工具箱的話,可以用Matlab里面的haltonset,faureset,sobolset函數(shù)實(shí)現(xiàn),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