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

蒲丰氏投针问题的模拟过程，随机数发生器也是自编的，以供大家参考和提出建议。谢谢。（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)

pi1=(1.0*n1)/(counter*1.0)*(2.0*L/a) write(*,*) 'this is PI:' write(*,*) pi

write(*,\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*x^2)在(0,1)间积分 format long;

syms x

a = sym(1/2);

f = exp(-a*x^2);

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); end

mju=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); end

mju=mju/n;

fprintf('Monte Carlo result:%1.18f.\\n',mju);

%=============生成HALTON序列========================

function result = halton( m,base,seeder )

%生成HALTON序列

% Check inputs

if nargin < 3

seeder = 0;

if nargin < 2

error('MATLAB:Halton:NotEnoughInputs',...

'Not enough input arguments. See Halton.'); end end

res=0;

n=length(base);

for i=1:m

for j=1:n

element=0;

temp=seeder+i;

k=1;

while temp>0

element(k)=rem(temp,base(j));

temp=fix(temp/base(j));

k=k+1;

end

res(i,j)= 0;

for k=1:length(element)

res(i,j)=res(i,j)+element(k)/(base(j)^k);

end

end

end

result=res;

本文来源：https://www.bwwdw.com/article/smbr.html