第7章 蒙特卡罗方法 附录.docx
《第7章 蒙特卡罗方法 附录.docx》由会员分享，可在线阅读，更多相关《第7章 蒙特卡罗方法 附录.docx（37页珍藏版）》请在冰豆网上搜索。
第7章蒙特卡罗方法附录
第7章附录
7.2.1均匀分布随机数
例题7.2.1计算程序
!
rand1.for
programrand1
implicitnone
realr
integern,c,x,i
open（5,file='rand1.txt'）
n=32768
c=889
x=13
doi=1,1000
x=c*x-n*int（c*x/n）
r=real（x）/（n-1）
write（5,'（f8.5）'）r
enddo
end
!
!
!
!
!
!
rand2.for!
!
!
!
!
programrand2
implicitnone
integer,parameter:
:
n=1000
integerix,i
realr
open（5,file='rand2.txt'）
ix=32765
doi=1,n
callrand（ix,r）
write（5,'（f8.6）'）r
enddo
endprogramrand2
subroutinerand（ix,r）
i=ix*259
ix=i-i/32768*32768
r=float（ix）/32768
return
end
7.2.3随机抽样
例题7.2.2计算程序
%例题7_2_2.m
figure
（1）;
set（gca,'FontSize',16）;
t=rand（1000,1）;
y=-log（t）;
z=exp（-y）;
plot（y,z,'.'）;
xlabel（'图7.2-2例题7.2.2-指数分布抽样'）
====================================================
例题7.2.5计算程序
!
例题7.2.5
programscores
parameter（nmax=10,mmax=13）
real（8）x（nmax）,y（nmax）,l（0:
nmax）,z（mmax）,ys（mmax）,r
integeri,j,k
datax/5.0,15.0,25.0,35.0,45.0,55.0,65.0,75.0,85.0,95.0/
datay/0,0,0,0,0.08,0.19,0.31,0.27,0.11,0.04/
open（2,file='scores_old.txt'）
open（5,file='scores_new.txt'）!
mmax个抽样学生成绩
open（7,file='scores_sample.txt'）
write（2,'（2f15.5）'）（x（i）,y（i）,i=1,nmax）
ix=32765
l（0）=0
doi=1,nmax
l（i）=l（i-1）+y（i）
enddo
doj=1,mmax
callrand（ix,r）
dok=1,nmax
if（r.le.l（k））goto11
enddo
11z（j）=x（k）
enddo
write（5,*）（z（i）,i=1,mmax）
ys=0
doi=1,mmax
k=z（i）/float（nmax）!
确定抽样学生所在的分数段
ys（k）=ys（k）+1.0!
统计每个分数段的学生数
enddo
doi=1,nmax
ys（i）=ys（i）/float（mmax）
write（7,*）x（i）,ys（i）
enddo
end
subroutinerand（ix,r）
integeri
real（8）r
i=ix*259
ix=i-i/32768*32768
r=float（ix）/32768
return
end
====================================================
7.3.1方程求根的M-C方法
例题7.3.1计算程序
!
例题3.1
programroots_MC1
implicitnone
realb,bmax,a,h,f,x
open（5,file='x.txt'）
b=0.0;bmax=4.0;h=0.3
4a=b
b=b+h
if（（f（a）*f（b））.gt.0.0）goto4
callsr（a,b,x）
write（5,'（f12.8）'）x
write（*,'（f12.8）'）x
if（b.lt.bmax）goto4
end
subroutinesr（a,b,x）
implicitnone
realr,xr,x,a,b,p,f
integerix,n,m,i
ix=32765;n=4000;m=0
doi=1,n
callrand（ix,r）
xr=a+（b-a）*r
if（f（xr）.le.0.0）m=m+1
enddo
p=real（m）/real（n）
if（f（a）.le.0.0）then
x=a+（b-a）*p
else
x=a+（b-a）*（1.0-p）
endif
return
end
functionf（x）
implicitnone
realx,f
!
f=exp（x）*（log（x）-x*x）
f=（（（x-10.0）*x+35.0）*x-50.0）*x+24.0
return
end
subroutinerand（ix,r）
i=ix*259
ix=i-i/32768*32768
r=float（ix）/32768
return
end
=========================================
例题7.3.2计算程序
!
!
!
!
!
!
!
!
!
!
!
!
roots_MC_2.f90!
!
!
!
!
!
!
!
!
!
!
programroots_MC_2
implicitnone
real（8）x,a,b,r,rp,rm,f,p
integerix,n,i
a=0.0;b=3.1415926/2.0
!
a=0.0;b=1.0
ix=32765
n=2000
rp=0.0;rm=1.0!
注意：
首先将最大赋值最小，将最小赋值最大
if（f（a）<0）then
doi=1,2000
callrand（ix,r）
if（f（a+（b-a）*r）<=0.0）then
callamax（r,rp）
else
callamin（r,rm）
endif
enddo
else
doi=1,2000
callrand（ix,r）
if（f（a+（b-a）*r）>=0.0）then
callamax（r,rp）
else
callamin（r,rm）
endif
enddo
endif
p=（rp+rm）/2.0
x=a+（b-a）*p
write（*,*）'x=',x
endprogramroots_MC_2
subroutineamax（x,xmax）
implicitnone
real（8）x,xmax
if（xmax.ge.x）then
return
else
xmax=x
endif
return
end
subroutineamin（x,xmin）
implicitnone
real（8）x,xmin
if（xmin.le.x）then
return
else
xmin=x
endif
return
end
functionf（x）
real（8）x
f=exp（-x*x*x）-tan（x）+800.0
!
f=x-exp（-x）
return
end
subroutinerand（ix,r）
real（8）r
integerix
i=ix*259
ix=i-i/32768*32768
r=float（ix）/32768
return
end
===========================================
7.3.2计算定积分的M-C方法
例题7.3.3计算程序
!
例题7.3.3
programint_MC1
real（8）f,s,a,b
a=2.5!
积分下限
b=8.4!
积分上限
callint_mc（a,b,s）
print*,'s=',s
end
functionf（x）
real（8）x
f=x*x+sin（x）
end
subroutineint_mc（a,b,s）
real（8）a,b,f,s,r,x
integerm,i
r=1.0
m=10000
s=0.0
doi=1,m
x=a+（b-a）*rand（r）
s=s+f（x）
enddo
s=s*（b-a）/real（m）
return
end
functionrand（r）
real（8）s,u,v,r
integerm
s=65536.0
u=2053.0
v=13849.0
m=r/s
r=r-m*s
r=u*r+v
m=r/s
r=r-m*s
rand=r/s
return
end
=================================================
例题7.3.4计算程序
!
例题3.4
programpi_MC1
integerm,n,s,ntotal
real（8）r,x,y,rand,iseed
open（5,file='pi.txt'）
iseed=1.0
n=100000
m=10000
s=0
dontotal=1,n
x=rand（iseed）
y=rand（iseed）
r=sqrt（x*x+y*y）
if（r<=1.0）s=s+1
if（mod（ntotal,m）.eq.0）then
print*,ntotal,4.0*float（s）/float（ntotal）
write（5,'（2x,i8,3x,f18.8）'）ntotal,4.0*float（s）/float（ntotal）
endif
enddo
end
functionrand（r）
real（8）s,u,v,r
integerm
s=65536.0
u=2053.0
v=13849.0
m=r/s
r=r-m*s
r=u*r+v
m=r/s
r=r-m*s
rand=r/s
return
end
%pi的计算（例题7.3.4）
clc;clear;formatlong;
i=1;count=0.;
n=5000;
x=rand（n,1）;
y=rand（n,1）;
figure
（1）;
set（gca,'FontSize',16）;
fori=1:
n
if（x（i）*x（i）+y（i）*y（i））<=1.0
plot（x（i）,y（i）,'r.'）;
count=count+1.;
else
holdon;
plot（x（i）,y（i）,'b.'）;
end
end
pi_val=4.*count/（i-1）
xlabel（'x'）;
ylabel（'y'）;
title（'图7.3-2计算\pi值示意图'）;
============================================================
例题7.3.5计算程序
!
例题3.5
programMint_MC1
implicitnone
real（8）c,xmax,rand
integerix,n,m,j,k
open（4,file='gama.txt'）
ix=32765
n=2000
dom=5,9
c=0.0
doj=1,n
xmax=rand（ix）
dok=2,m
xmax=max（xmax,rand（ix））
enddo
c=c+dlog（dabs（m*dlog（xmax）））
enddo
c=c/n
write（4,100）m,c
write（*,100）m,c
enddo
100format（2x,'gama（',i2,'）=',f18.8）
end
realfunctionrand（ix）
integeri,ix
i=ix*259
ix=i-i/32768*32768
rand=real（ix）/32768.0
return
end
=========================================
例题7.3.6计算程序
!
!
!
!
int_MC2.f90!
!
!
!
!
!
!
programint_MC2
integerm,n,s,ntotal
real（8）r,x,y,rand,iseed
open（5,file='pi.txt'）
iseed=1.0
n=100000
m=10000
s=0
dontotal=1,n
x=rand（iseed）
y=rand（iseed）
r=sqrt（x*x+y*y）
if（r<=1.0）s=s+1
if（mod（ntotal,m）.eq.0）then
print*,ntotal,4.0*float（s）/float（ntotal）
write（5,'（2x,i8,3x,f18.8）'）ntotal,4.0*float（s）/float（ntotal）
endif
enddo
end
functionrand（r）
real（8）s,u,v,r
integerm
s=65536.0
u=2053.0
v=13849.0
m=r/s
r=r-m*s
r=u*r+v
m=r/s
r=r-m*s
rand=r/s
return
end
/*Mint_MC2.c例题3.6*/
#include
#definef（x,y）x*y*y
staticdoublexn=1;
doublerand（）
{doublelambda=16807,p=2147483647;
xn=fmod（lambda*xn,p）;
returnxn/p;}
voidsrand（doubleseed）
{xn=seed;}
main（）
{inti,n=15000;doublex,y,s=0;
srand（4）;
for（i=0;i{x=rand（）+2/3.0;y=（3*x-2）*rand（）+（2-x）;
s=s+（3*x-2）*f（x,y）;}
s=s/n;
printf（"N=%8dV=%10.6f\n",n,s）;
}
例题7.3.7计算程序
/*Mint_MC3.c例题.3.7*/
#include
staticdoublexn=1;
doublerand（）
{doublelambda=16807,p=2147483647;
xn=fmod（lambda*xn,p）;
returnxn/p;}
voidsrand（doubleseed）
{xn=seed;}
main（）
{inti,n=10000,m=0;doublex,y,z,v=0;
srand（256）;
for（i=0;i{x=-0.5+rand（）;y=-0.5+rand（）;z=-0.5+rand（）;
if（（sqrt（x*x+y*y+z*z）<=0.5）&&（sqrt（x*x+y*y）>=0.3））
m=m+1;}
v=1.*（double）m/n;
printf（"N=%8dV=%10.6f\n",n,v）;
}
======================================================
7.3.3MC方法求解Laplace方程
例题7.3.8计算程序
!
programWalk_MC1.for
!
蒙特卡罗方法求解2d_Laplaciandifferenceequation
parameter（imax=50,jmax=50,pi=3.1415926,kmax=10000）
dimensionT（imax,jmax）,a（imax,jmax）,
&qx（imax,jmax）,qy（imax,jmax）,q（imax,jmax）,an（imax,jmax）
reallamda,Tu,Td,Tl,Tr,kpx,kpy,dx,dy,r
open（2,file='walk2.dat'）
open（3,file='q2.dat'）
lamda=1.5;Tu=100.;Tl=75.;Tr=50.;Td=0.;dx=0.1;dy=0.1;ix=32765
kpx=-0.49/2./dx;kpy=-0.49/2./dy
T（:
:
）=0.0
T（1,1）=0.5*（Td+Tl）;T（1,jmax）=0.5*（Tl+Tu）
T（imax,jmax）=0.5*（Tu+Tr）;T（imax,1）=0.5*（Tr+Td）
T（1,2:
jmax-1）=Tl;T（imax,2:
jmax-1）=Tr
T（2:
imax-1,1）=Td;T（2:
imax-1,jmax）=Tu
doi=2,imax-1
print*,'i=',i
ii=i
doj=2,jmax-1
jj=j
do20k=1,kmax
a（i,j）=T（i,j）
6callrand（ix,r）
if（r.le.0.25）then
ii=ii+1
else
if（r.le.0.5）then
jj=jj-1
else
if（r.le.0.75）then
ii=ii-1
else
jj=jj+1
endif
endif
endif
if（ii.eq.1.or.ii.eq.imax）goto18
if（jj.eq.1.or.jj.eq.jmax）goto19
goto6
18T（i,j）=T（ii,j）
goto17
19T（i,j）=T（i,jj）
17T（i,j）=a（i,j）+T（i,j）
20continue
T（i,j）=T（i,j）/float（kmax）
enddo
enddo
doj=1,jmax
write（2,222）（T（i,j）,i=1,imax）
enddo
222format（1x,50（2x,f14.10））
doi=2,imax-1
doj=2,jmax-1
qx（i,j）=kpx*（T（i+1,j）-T（i-1,j））
qy（i,j）=kpy*（T（i,j+1）-T（i,j-1））
q（i,j）=sqrt（qx（i,j）*qx（i,j）+qy（i,j）*qy（i,j））
if（qx（i,j）.gt.0.）then
an（i,j）=atan（qy（i,j）/qx（i,j））
else
an（i,j）=atan（qy（i,j）/qx（i,j））+pi
endif
enddo
enddo
doj=2,jmax-1
doi=2,imax-1
!
write（3,333）dx*（i-1）,dy*（j-1）,an（i,j）,q（i,j）
write（3,333）dx*（i-1）,dy*（j-1）,qx（i,j）,qy（i,j）
enddo
enddo
333format（1x,4（2x,f18.12））
end
subroutinerand（ix,r）
i=ix*259
ix=i-i/32768*32768
r=float（ix）/32768
return
end
======================================================
7.3.4核链式反应的模拟
!
!
!
!
!
!
MC_1.f90!
!
!
!
!
programmc_1
implicitnone
integeri,j,k,nmax,kmax
real（8）x,n,m,dm,t,pi,a,b,nin,rand,f,xi（50）,fi（50）
real（8）r（9）,x0,y0,z0,x1,y1,z1,p,c,s,d
open（2,file='k_m.txt'）
x=1.0;n=10000.0;m=0.98
dm=0.02;t=1.0;pi=3.14159265
nmax=n;kmax=50
dok=1,kmax
m=m+dm
xi（k）=m
a=（m*t）**（1.0/3.0）
b=（m/（t*t））**（1.0/3.0）
nin=0.0
doj=1,nmax
doi=1,9
r（i）=rand（x）
enddo
x0=a*（r
（1）-0.5）
y0=a*（r
（2）-0.5）
z0=b*（r（3）-0.5）
do1i=1,2
p=2.0*pi*r（2*i+2）
c=2.0*（r（2*i+3）-0.5）
s=sqrt（1.0-c*c）
d=r（i+7）
x1=x0+d*s*cos（p）
y1=y0+d*s*sin（p）
z1=z0+d*c
if（abs（x1）.gt.0.5*a）goto1
if（abs（y1）.gt.0.5*a）goto1
if（abs（z1）.gt.0.5*b）goto1
nin=nin+1.0
1continue
enddo
f=nin/n
fi（k）=f
write（2,'（5x,d13.8,5x,d13.8）'）m,f
enddo
endprogrammc_1
rea
升级会员 