《第7章-蒙特卡罗方法-(附录)(共25页).doc》由会员分享，可在线阅读，更多相关《第7章-蒙特卡罗方法-(附录)(共25页).doc（25页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精选优质文档-倾情为你奉上第7章 附录7.2.1 均匀分布随机数例题7.2.1计算程序! rand1.forprogram rand1implicit nonereal rinteger n,c,x,i open(5,file=rand1.txt)n = 32768c = 889x = 13 do i = 1,1000 x = c*x-n*int(c*x/n) r = real(x)/(n-1) write(5,(f8.5) rend doend!rand2.for!program rand2implicit noneinteger, parameter : n=1000integer ix,
2、ireal ropen(5,file=rand2.txt)ix=32765do i=1,n call rand(ix,r) write(5,(f8.6) rend doend program rand2subroutine rand(ix,r)i=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend7.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.
3、2-2 例题7.2.2-指数分布抽样)=例题7.2.5计算程序! 例题7.2.5 program scoresparameter(nmax=10,mmax=13)real(8) x(nmax),y(nmax),l(0:nmax),z(mmax),ys(mmax),rinteger i,j,kdata x/5.0,15.0,25.0,35.0,45.0,55.0,65.0,75.0,85.0,95.0/data y/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
4、) ! mmax个抽样学生成绩open(7,file=scores_sample.txt) write(2,(2f15.5) (x(i),y(i),i=1,nmax)ix=32765 l(0)=0do i=1,nmax l(i)=l(i-1)+y(i)end do do j=1,mmax call rand(ix,r) do k=1,nmax if(r.le.l(k) goto 11 end do11 z(j)=x(k) end do write(5,*) (z(i),i=1,mmax)ys=0do i=1,mmax k=z(i)/float(nmax) ! 确定抽样学生所在的分数段 ys(k
5、)=ys(k)+1.0 ! 统计每个分数段的学生数end do do i=1,nmax ys(i)=ys(i)/float(mmax) write(7,*) x(i),ys(i)end doendsubroutine rand(ix,r)integer ireal(8) ri=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend=7.3.1 方程求根的M-C方法例题7.3.1计算程序! 例题3.1program roots_MC1 implicit nonereal b,bmax,a,h,f,xopen(5,file=x.txt) b=0.0;
6、bmax=4.0; h=0.34 a=b b=b+hif(f(a)*f(b).gt.0.0) goto 4 call sr(a,b,x)write(5,(f12.8) x write(*,(f12.8) x if(b.lt.bmax) goto 4end subroutine sr(a,b,x)implicit nonereal r,xr,x,a,b,p,finteger ix,n,m,iix=32765; n = 4000; m = 0 do i =1,n call rand(ix,r) xr=a+(b-a)*r if(f(xr).le.0.0) m = m+1end dop = real(
7、m)/real(n) if(f(a).le.0.0) then x=a+(b-a)*pelse x=a+(b-a)*(1.0-p)end ifreturnendfunction f(x) implicit nonereal x,f!f=exp(x)*(log(x)-x*x)f=(x-10.0)*x+35.0)*x-50.0)*x+24.0returnendsubroutine rand(ix,r)i=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend=例题7.3.2计算程序! roots_MC_2.f90 !program roots_MC_2
8、implicit nonereal(8) x,a,b,r,rp,rm,f,pinteger ix,n,i a=0.0; b=3./2.0! a=0.0;b=1.0ix=32765n=2000 rp=0.0; rm=1.0 !注意：首先将最大赋值最小，将最小赋值最大if(f(a)0) then do i=1,2000 call rand(ix,r) if(f(a+(b-a)*r)=0.0) then call amax(r,rp) else call amin(r,rm) end if end do end if p=(rp+rm)/2.0 x=a+(b-a)*p write(*,*) x=,x
9、end program roots_MC_2subroutine amax(x,xmax)implicit nonereal(8) x,xmax if(xmax.ge.x) then return else xmax=x end ifreturnendsubroutine amin(x,xmin)implicit nonereal(8) x,xmin if(xmin.le.x) then return else xmin=x end ifreturnendfunction f(x)real(8) xf=exp(-x*x*x)-tan(x)+800.0! f=x-exp(-x)returnend
10、subroutine rand(ix,r)real(8) rinteger ixi=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend=7.3.2 计算定积分的M-C方法例题7.3.3计算程序 ! 例题7.3.3program int_MC1real(8) f,s,a,ba = 2.5 ! 积分下限b = 8.4 ! 积分上限call int_mc(a,b,s)print *,s = , send function f(x)real(8) xf=x*x+sin(x)endsubroutine int_mc(a,b,s)real(8) a,b,f
11、,s,r,xinteger m,ir = 1.0m = 10000s = 0.0do i = 1,m x = a+(b-a)*rand(r) s = s+f(x)end dos=s*(b-a)/real(m)returnendfunction rand(r)real(8) s,u,v,rinteger ms=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.4 program pi_MC1integer m,n,s,ntotalreal(8) r,x,y,rand,
12、iseedopen(5,file=pi.txt)iseed = 1.0n = m = 10000s = 0do ntotal=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) endifend doend function rand(r)real(8) s,u,
13、v,rinteger ms=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; format long;i=1;count=0.;n = 5000;x = rand(n,1);y = rand(n,1);figure(1);set(gca,FontSize,16);for i=1:n if (x(i)*x(i)+y(i)*y(i)=1.0 plot(x(i),y(i),r.); count=count+1.; else hold on; plot(x
14、(i),y(i),b.); endendpi_val = 4.*count/(i-1)xlabel(x);ylabel(y);title(图7.3-2 计算pi 值示意图);=例题7.3.5计算程序! 例题3.5program Mint_MC1implicit nonereal(8) c,xmax,randinteger ix,n,m,j,kopen(4,file=gama.txt)ix=32765n=2000do m=5,9 c=0.0 do j=1,n xmax = rand(ix) do k=2,m xmax = max(xmax,rand(ix) end do c=c+dlog(dab
15、s(m*dlog(xmax) end do c=c/n write(4,100) m,c write(*,100) m,cend do100 format(2x,gama(,i2,)=,f18.8) end real function rand(ix) integer i,ixi=ix*259ix=i-i/32768*32768rand = real(ix)/32768.0returnend=例题7.3.6计算程序 ! int_MC2.f90 ! program int_MC2integer m,n,s,ntotalreal(8) r,x,y,rand,iseedopen(5,file=pi.
16、txt)iseed = 1.0n = m = 10000s = 0do ntotal=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) endifend doend function rand(r)real(8) s,u,v,rinteger ms=65536.
17、0u=2053.0v=13849.0m=r/sr=r-m*sr=u*r+vm=r/sr=r-m*srand=r/sreturnend/* Mint_MC2.c 例题3.6 */#include#define f(x,y) x*y*ystatic double xn=1;double rand() double lambda=16807, p=;xn=fmod(lambda*xn,p);return xn/p; void srand(double seed) xn=seed; main() int i,n=15000; double x, y, s=0;srand(4);for(i=0; in;
18、 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=%8d V=% 10.6fn, n, s);例题7.3.7计算程序/* Mint_MC3.c 例题.3.7 */#includestatic double xn=1;double rand() double lambda=16807, p=;xn=fmod(lambda*xn,p);return xn/p; void srand(double seed) xn=seed; main() int i,n=10000,m=0; double
19、x, y, z, v=0;srand(256);for(i=0; in; i+) x=-0.5+rand(); y=-0.5+rand(); z = -0.5+rand();if(sqrt(x*x+y*y+z*z)=0.3)m = m+1;v = 1.*(double)m/n;printf(N=%8d V=% 10.6fn, n, v);=7.3.3 MC方法求解Laplace方程例题7.3.8计算程序! program Walk_MC1.for! 蒙特卡罗方法求解2d_Laplacian difference equationparameter(imax=50,jmax=50,pi=3.,k
20、max=10000)dimension T(imax,jmax),a(imax,jmax), & qx(imax,jmax),qy(imax,jmax),q(imax,jmax),an(imax,jmax)real lamda,Tu,Td,Tl,Tr,kpx,kpy,dx,dy,ropen(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=32765kpx=-0.49/2./dx;kpy=-0.49/2./dyT(:,:)=0.0T(1,1)=0.5*(Td+Tl
21、) ;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)=Tudo i=2,imax-1 print *,i=,i ii=i do j=2,jmax-1 jj=j do 20 k=1,kmax a(i,j)=T(i,j)6 call rand(ix,r) if(r.le.0.25) then ii=ii+1 else if(r.le.0.5) then jj=jj-1
22、else if(r.le.0.75) then ii=ii-1 else jj=jj+1 end if end if end if if(ii.eq.1.or.ii.eq.imax) goto 18 if(jj.eq.1.or.jj.eq.jmax) goto 19 goto 618 T(i,j)=T(ii,j) goto 1719 T(i,j)=T(i,jj)17 T(i,j)=a(i,j)+T(i,j)20 continue T(i,j)=T(i,j)/float(kmax) end doend dodo j=1,jmax write(2,222) (T(i,j),i=1,imax)end
23、 do222format(1x,50(2x,f14.10)do i=2,imax-1 do j=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 end if end do end dodo j=2,jmax-1 do i=2,imax-1!
24、 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) end doend do333 format(1x,4(2x,f18.12)end subroutine rand(ix,r)i=ix*259ix=i-i/32768*32768r=float(ix)/32768returnend=7.3.4 核链式反应的模拟!MC_1.f90!program mc_1implicit noneinteger i,j,k,nmax,kmaxreal(8) x,n,m,dm,t
25、,pi,a,b,nin,rand,f,xi(50),fi(50)real(8) r(9),x0,y0,z0,x1,y1,z1,p,c,s,dopen(2,file=k_m.txt)x = 1.0 ;n = 10000.0; m = 0.98dm= 0.02; t = 1.0; pi= 3.nmax=n; kmax=50do k=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 do j=1,nmax do i=1,9 r(i)=rand(x) end do x0=a*(r(1)-0.5) y0=a*
26、(r(2)-0.5) z0=b*(r(3)-0.5) do 1 i=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) goto 1 if(abs(y1).gt.0.5*a) goto 1 if(abs(z1).gt.0.5*b) goto 1 nin=nin+1.01 continue end do f=nin/n fi(k)=f write(2,(5x,d13.8,5x,d13.
27、8) m,fend doend program mc_1real function rand(x)real(8) s,u,v,xs=65536.0u=2053.0v=13849.0m=x/sx=x-m*sx=u*x+vm=x/sx=x-m*srand=x/sreturnend% piot_k_m %z=load(k_m.txt);x=z(:,1);y=z(:,2);plot(x,y,ro);grid on;xlabel(M);ylabel(k);title(裂变倍增系数与质量之间关系);= 7.3.5 关于中子贯穿概率问题/* MC1.c */#includestatic double xn=
28、1;double rand() double lambda=16807,p=;xn=fmod(lambda*xn,p);return xn/p; void srand(double seed) xn=seed; main() int i,n= 10000,k,na=0,np=0,nr=0; double pi=3.,s,theta,pa,pp,pr; srand(256); for(i=1;i=n;i+) k=0;s=0; do theta=2*pi*rand(); s+=cos(theta);k+; if(s5)np+;break; if(k10)na+;break; while(1); p
29、p=(double)np/n* 100; pa=(double)na/n* 100; pr=(double)nr/n* 100; printf(PP=%6.2f%PA=%6.2f%PR=%6.2f%n,pp,pa,pr);= 例题7.3.9计算程序! MC_2.f90 !program real(8) h(3) ix=32765; r1=15.0; r2=16.0ab=8.33; aniu=0.1; b=0.4data h/25.0,115.0,1015.0/ do j=1,3 phi0 = acos(sqrt(h(j)*h(j)-r2*r2)/h(j) d = 0.0; phi = 0.0
30、do n = 1,10000 call rand(ix,x); call rand(ix,y) s=sqrt(2.0*x-x*x); c=1.0-x fi=phi0*y; oa=h(j)*sin(fi) ap=h(j)*cos(fi) ; ar2=sqrt(r2*r2-oa*oa) if(oa.ge.r1) then al=2.0*ar2/c else ar1=sqrt(r1*r1-oa*oa); al0=2.0*(ar2-ar1)/c al1g=int(ap+ar1)*s/c+0.5*ab)/ab) al1s=int(ap-ar1)*s/c+0.5*ab)/ab) al=al0+b/s*(a
31、l1g-al1s) end if dphi=phi0/3.1416/aniu*(1.0-exp(-aniu*al) phi=phi+dphi end do phi=phi/float(n) write(*,(2f14.7) h(j),phi end do endSUBROUTINE RAND(IX,R)I=IX*259IX=I-I/32768*32768R=FLOAT(IX)/32768RETURNEND= 7.3.6 其他例子1. 随机行走问题!RW.f90!PROGRAM Random_Walk!-! Monte Carlo Simulation of 1D Random Walk! Th
32、e walker begins at the origin (x=0) and first! step is choosen at random to be either to the! right or left each with probability 1/2 as in Figure.! -a 0 +a! x! ! The program calculates Average Number of Steps! for the walker to be outside of the region -a,+a.! Jan 2004! Ahmet Bingul, !-IMPLICIT NON
33、EINTEGER : I,J,N,X,Step,Sum,A=3REAL : R,AvrS! initiate random number generator CALL RANDOM_SEED DO I=1,6 N = 10*I ! Number of Experiment Sum = 0 ! Sum of steps DO J=1,N Step = 0 ! step counter X = 0 ! current location of the walker !- Do an experiment -! DO ! ! CALL RANDOM_NUMBER(R) ! Generate a random number ! IF(R A ) THEN ! Is the walker outside? Sum = Sum + Step ! Sum the steps EXIT ! END IF ! ! END DO ! !-! END DO AvrS = REAL(Sum)/N ! Average number of steps PRINT (I15,F10.6),N,AvrS ! output N vs Avrs END DOEND PROGRAM Random_Walk