第7章蒙特卡罗方法(附录)

1第7章附录7.2.1均匀分布随机数例题7.2.1计算程序!rand1.forprogramrand1implicitnonerealrintegern,c,x,iopen(5,file='rand1.txt')n=32768c=889x=13doi=1,1000x=c*x-n*int(c*x/n)r=real(x)/(n-1)write(5,'(f8.5)')renddoend!!!!!!rand2.for!!!!!programrand2implicitnoneinteger,parameter::n=1000integerix,irealropen(5,file='rand2.txt')ix=32765doi=1,ncallrand(ix,r)write(5,'(f8.6)')renddoendprogramrand2subroutinerand(ix,r)i=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend27.2.3随机抽样例题7.2.2计算程序%例题7_2_2.mfigure(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.5programscoresparameter(nmax=10,mmax=13)real(8)x(nmax),y(nmax),l(0:nmax),z(mmax),ys(mmax),rintegeri,j,kdatax/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=32765l(0)=0doi=1,nmaxl(i)=l(i-1)+y(i)enddodoj=1,mmaxcallrand(ix,r)dok=1,nmaxif(r.le.l(k))goto11enddo11z(j)=x(k)enddowrite(5,*)(z(i),i=1,mmax)ys=0doi=1,mmaxk=z(i)/float(nmax)!确定抽样学生所在的分数段3ys(k)=ys(k)+1.0!统计每个分数段的学生数enddodoi=1,nmaxys(i)=ys(i)/float(mmax)write(7,*)x(i),ys(i)enddoendsubroutinerand(ix,r)integerireal(8)ri=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend====================================================7.3.1方程求根的M-C方法例题7.3.1计算程序!例题3.1programroots_MC1implicitnonerealb,bmax,a,h,f,xopen(5,file='x.txt')b=0.0;bmax=4.0;h=0.34a=bb=b+hif((f(a)*f(b)).gt.0.0)goto4callsr(a,b,x)write(5,'(f12.8)')xwrite(*,'(f12.8)')xif(b.lt.bmax)goto4endsubroutinesr(a,b,x)implicitnonerealr,xr,x,a,b,p,fintegerix,n,m,iix=32765;n=4000;m=0doi=1,ncallrand(ix,r)4xr=a+(b-a)*rif(f(xr).le.0.0)m=m+1enddop=real(m)/real(n)if(f(a).le.0.0)thenx=a+(b-a)*pelsex=a+(b-a)*(1.0-p)endifreturnendfunctionf(x)implicitnonerealx,f!f=exp(x)*(log(x)-x*x)f=(((x-10.0)*x+35.0)*x-50.0)*x+24.0returnendsubroutinerand(ix,r)i=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend=========================================例题7.3.2计算程序!!!!!!!!!!!!roots_MC_2.f90!!!!!!!!!!!programroots_MC_2implicitnonereal(8)x,a,b,r,rp,rm,f,pintegerix,n,ia=0.0;b=3.1415926/2.0!a=0.0;b=1.0ix=32765n=2000rp=0.0;rm=1.0!注意：首先将最大赋值最小，将最小赋值最大if(f(a)0)thendoi=1,2000callrand(ix,r)if(f(a+(b-a)*r)=0.0)thencallamax(r,rp)5elsecallamin(r,rm)endifenddoelsedoi=1,2000callrand(ix,r)if(f(a+(b-a)*r)=0.0)thencallamax(r,rp)elsecallamin(r,rm)endifenddoendifp=(rp+rm)/2.0x=a+(b-a)*pwrite(*,*)'x=',xendprogramroots_MC_2subroutineamax(x,xmax)implicitnonereal(8)x,xmaxif(xmax.ge.x)thenreturnelsexmax=xendifreturnendsubroutineamin(x,xmin)implicitnonereal(8)x,xminif(xmin.le.x)thenreturnelsexmin=xendifreturnendfunctionf(x)real(8)x6f=exp(-x*x*x)-tan(x)+800.0!f=x-exp(-x)returnendsubroutinerand(ix,r)real(8)rintegerixi=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend===========================================7.3.2计算定积分的M-C方法例题7.3.3计算程序!例题7.3.3programint_MC1real(8)f,s,a,ba=2.5!积分下限b=8.4!积分上限callint_mc(a,b,s)print*,'s=',sendfunctionf(x)real(8)xf=x*x+sin(x)endsubroutineint_mc(a,b,s)real(8)a,b,f,s,r,xintegerm,ir=1.0m=10000s=0.0doi=1,mx=a+(b-a)*rand(r)s=s+f(x)enddos=s*(b-a)/real(m)7returnendfunctionrand(r)real(8)s,u,v,rintegerms=65536.0u=2053.0v=13849.0m=r/sr=r-m*sr=u*r+vm=r/sr=r-m*srand=r/sreturnend=================================================例题7.3.4计算程序!例题3.4programpi_MC1integerm,n,s,ntotalreal(8)r,x,y,rand,iseedopen(5,file='pi.txt')iseed=1.0n=100000m=10000s=0dontotal=1,nx=rand(iseed)y=rand(iseed)r=sqrt(x*x+y*y)if(r=1.0)s=s+1if(mod(ntotal,m).eq.0)thenprint*,ntotal,4.0*float(s)/float(ntotal)write(5,'(2x,i8,3x,f18.8)')ntotal,4.0*float(s)/float(ntotal)endifenddoendfunctionrand(r)real(8)s,u,v,r8integerms=65536.0u=2053.0v=13849.0m=r/sr=r-m*sr=u*r+vm=r/sr=r-m*srand=r/sreturnend%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:nif(x(i)*x(i)+y(i)*y(i))=1.0plot(x(i),y(i),'r.');count=count+1.;elseholdon;plot(x(i),y(i),'b.');endendpi_val=4.*count/(i-1)xlabel('x');ylabel('y');title('图7.3-2计算\pi值示意图');============================================================例题7.3.5计算程序!例题3.5programMint_MC1implicitnonereal(8)c,xmax,randintegerix,n,m,j,k9open(4,file='gama.txt')ix=32765n=2000dom=5,9c=0.0doj=1,nxmax=rand(ix)dok=2,mxmax=max(xmax,rand(ix))enddoc=c+dlog(dabs(m*dlog(xmax)))enddoc=c/nwrite(4,100)m,cwrite(*,100)m,cenddo100format(2x,'gama(',i2,')=',f18.8)endrealfunctionrand(ix)integeri,ixi=ix*259ix=i-i/32768*32768rand=real(ix)/32768.0returnend=========================================例题7.3.6计算程序!!!!int_MC2.f90!!!!!!!programint_MC2integerm,n,s,ntotalreal(8)r,x,y,rand,iseedopen(5,file='pi.txt')iseed=1.0n=100000m=10000s=0dontotal=1,nx=rand(iseed)y=rand(iseed)r=sqrt(x*x+y*y)if(r=1.0)s=s+110if(mod(ntotal,m).eq.0)thenprint*,ntotal,4.0*float(s)/float(ntotal)write(5,'(2x,i8,3x,f18.8)')ntotal,4.0*float(s)/float(ntotal)endifenddoendfunctionrand(r)real(8)s,u,v,rintegerms=65536.0u=2053.0v=13849.0m=r/s