计算物理实验考试及答案.doc

07级计算物理实验期末考试试题解答姓名：马瑞松 学号：020112007016 班级：物理学2007级 成绩： 1、 用fortran 语言编程采用两种方法（高等数学方法和montecarlo方法）计算。答： 1、高等数学方法：program main parameter(pi=3.1415926) Double precision x,f1,f2,s,si,h f(x)=(sin(x)*4.0/x*3.0 h=2*pi/20110.0 s=0.0 f1=0.0 do 10 i=1,20110 f2=f(i*h) si=(f1+f2)*h/2.0 s=s+si f1=f210 continue print 100,s100 format (1x,f8.5) End 此程序经调试无误，输出结果为0.68860，以jifen.f保存在所提交所提交压缩包里。 2、蒙特卡洛方法： program main Double precision x,f,total,mean,result parameter(pi=3.1415926) f(x)=(sin(x)*4.0/x*3.0 total=0.0 do 10 i=1,2010000 x=2*pi*rand() total=f(x)+total10 continue mean=total/2010000 result=mean*2*pi print 100,result100 format (1x,result=,f7.5) End此程序输出结果为0.68860，以mtkljifen.f保存在所提交所提交压缩包里。2、 用Gear predictor-corrector算法数值求解正方形区域 的拉普拉斯方程答：所编程序如下所示：subroutine posor(n,m,n1,m1,u0,u1,u,h1,h2,f,l,maxl,int,w0,eps)dimension u0(m,n),u1(m,n),u(m,n),f(m1,n1)if(int.gt.0) goto 11l=0do 1 i=2,m1 do 1 j=2,n1 u0(i,j)=0.0 u1(i,j)=0.01 continue do 2 i=1,m u(i,1)=u0(i,1) u(i,n)=u0(i,n)2 continuedo 3 j=1,n u(1,j)=u0(1,j) u(m,j)=u0(m,j)3continue4 do 5 i=2,m1 do 5 j=2,n1 u(i,j)=(u0(i+1,j)+u(i-1,j)*h1*2+ $ (u0(i,j+1)+u(i,j-1)*h2*2-f(i,j)* $ (h1*h2)*2)/(2.0*(h1*2+h2*2)5 continue if(l.le.2) goto 8 p1=0.0p2=0.0do 6 i=1,mdo 6 j=1,np1=p1+abs(u(i,j)-u0(i,j)p2=p2+abs(u0(i,j)-u1(i,j)6 continuew1=p1/p2if(l.le.3) goto 7if(w1.gt.1.0) pause eror resultif(abs(w1-w0).lt.eps) goto 107 wo=w18l=l+1do 9 i=1,mdo 9 j=1,nu1(i,j)=u0(i,j)u0(i,j)=u0(i,j)9continuegoto 410w0=2.0/(1.0+sqrt(1.0-w1*w1)goto 1511l=016do 12 i=2,m1do 12 j=2,n1u(i,j)=u0(i,j)*(1.0-w0)+(u0(i+1,j)+ $ u(i-1,j)*h1*2+(u0(i,j+1)+u(i,j-1) * *h2*2-f(i,j)*(h1*h2)*2)*w0/(2.0*(h1*2+h2*2)12continueeror=0.0do 14 i=1,mdo 14 j=1,neror=eror+abs(u0(i,j)-u(i,j)14u0(i,j)=u(i,j)l=l+1if(l.gt.maxl) goto 15if(eror.gt.eps) goto 1615returnendprogram tposordimension u0(17,17),u1(17,17),u(17,17),f(16,16)eps=1.e-5n=17m=17n1=n-1m1=m-1h1=1.0/n1h2=1.0/m1do 1 i=1,mu0(i,1)=0.0u0(i,1)=0.01continuedo 2 j=2,n1u0(1,j)=1.0u0(m,j)=1.02continuedo 3 i=2,m1do 3 j=2,n1 f(i,j)=0.03continuecall posor(m,n,n1,m1,u0,u1,u,h1,h2,f,l,1000,0,w0,eps)write(*,30) l,w030format(2x,l=,i5,3x,w0=,f15.7)call posor(m,n,n1,m1,u0,u1,u,h1,h2,f,l,1000,1,w0,eps)write(*,10) l10format(2x,l=,i5)do 4 j=2,n1write(*,20) j,(u(i,j),i=2,m1/2+1)20format(1x,i2,8f5.3)4continuestopEnd此程序经过调试后以laplace.f保存在所提交所提交压缩包里。3、 十个氢原子在300K温度下其初速度满足麦克斯韦速率分布，即正比于，用fortran语言编程得到其初速度。答：所编程序如下所示：program main Double precision FF,v,zmass,t,k,f dimension x(10) parameter(k=1.3806505D-23,pi=3.1415926) FF(v,zmass,t)=(zmass/(2*pi*k*t)*(3.0/2) $ *exp(-zmass*(v*2)/(2*k*t) print*,please input the mass and temperture read*,zmass,t zmass=zmass*1.6726231D-27 v=0.0 fmax=FF(sqrt(2*k*t/zmass),zmass,t) cc=sqrt(2*k*t/zmass) print*,fmax=,fmax, vp= ,cc i=1 do 10 while(i.le.10) v=rand(0)*5000 f=rand(0)*(fmax) if(f.le.FF(v,zmass,t)then x(i)=v i=i+1 endif 10 continue print*,(x(i),i=1,10) End由键盘输入1 300，可以得到10个符合麦克斯韦分布的速度，此程序以maikesiwei.f保存在所提交压缩包内。4、 十个氢原子在模拟盒中运动，如果盒中有1000米/秒得风自左向右吹，用fortran语言编程计算左右器壁的压力差，适当选取模拟盒大小和温度，原子与器壁的碰撞是完全弹性的。答：所编程序如下所示： program main double precision ff,v,zmass,t,k,f,p1,p2 dimension x(10) parameter(k=1.3806505d-23,pi=3.1415926) ff(v,zmass,t)=(zmass/(2*pi*k*t)*(3.0/2) $ *exp(-zmass*(v*2)/(2*k*t) print*,please input the mass and temperture read*,zmass,t zmass=zmass*1.6726231d-27 v=0.0 fmax=ff(sqrt(2*k*t/zmass),zmass,t) i=1 do 10 while(i.le.10) v=rand(0)*5000 f=rand(0)*(fmax) if(f.le.ff(v,zmass,t)then x(i)=v i=i+1 endif 10 continue print*,(x(i),i=1,10) p1=0.0 do 20 i=1,10 p1=1.0/1D-21*zmass*(x(i)+1000)*2.0+p1 20 continue p2=0.0 do 30 i=1,10 p2=1.0/1d-21*zmass*(x(i)-1000)*2.0+p2 30 continue print *,(p1-p2),Pa print 100,(p1-p2)*1D-14100 format (1x,deltaF=,e8.3,N) end此程序令模拟盒为长度是1D-7m的正方体，由键盘输入氢原子相对原子质量1以及自行设定的温度，其初速度符合麦克斯韦速率分布，输出结果为左右器壁的压强差与压力差。调试后程序以press.f保存在所提交压缩包内。五、编写一个三维，元胞尺寸为L3的周期边界条件计算程序。答：所编程序如下所示：program main real pos(10000,3),v(1,3)parameter(boxl=1)gap=5E-6n=0 data (pos(1,j),j=1,3)/0,0,0/do 10,j=1,3 v(1,j)=100.0*rand(0)-5010continuedo 100,i=2,10000 pos(i,1)=pos(i-1,1)+v(1,1)*gap pos(i,2)=pos(i-1,2)+v(1,2)*gappos(i,3)=pos(i-1,3)+v(1,3)*gapcall bound(pos,i,n)100continueprint*,Please input which step you want to know!read*,iprint*,x= y= z= print*,(pos(i,j),j=1,3)print*,The total use of boundary method,nEndsubroutine bound(pos,i,n)real pos(10000,3)parameter(boxl=1)if (pos(i,1).lt.-boxl) thenpos(i,1)=pos(i,1)+boxl*2in=in+1else if(pos(i,1).gt.boxl) thenpos(i,1)=pos(i,1)-boxl*2in=in+1endifif(pos(i,2).lt.-boxl) thenpos(i,2)=pos(i,2)+boxl*2in=in+1else if(pos(i,2).gt.boxl) thenpos(i,2)=pos(i,2)-boxl*2in=in+1endifif(pos(i,3).lt.-boxl) thenpos(i,3)=pos(i,3)+boxl*2in=in+1else if(pos(i,3).gt.boxl) thenpos(i,3)=pos(i,3)-boxl*2in=in+1endifif(in.gt.0) n=n+1End此程序假设元胞长度L为1，模拟一个原子在元胞内采用周期性边界条件的运动。原子的初始位置以及初始速度随机给出，由键盘输入欲想知道的步数，即可得到该原子此时的位置以及周期性边界条件采用的次数。此程序以boundary.f保存在所提交压缩包内。6、 试做总能量固定的单原子系统的分子动力学模拟。元胞为，划分为的正方形网格。元胞内原子数。原子质量。位势为Lenard-Jones势，其中，边界条件为周期性边界条件，初始位置是随机分布在正则节点上，初始速度为按-1，1随机分布。分子动力学模拟步长取为，模拟100-200步后原子的速度分布和位置分布如何？答：所编程序如下所示：program main dimension pos(64,3,200),v(64,3,200),a(64,3,200)data a/38400*0.0/call begin(pos,v,a)call getac(pos,a,1)do nstep=2,200call getr(pos,v,a,nstep)call bound(pos,nstep)call getac(pos,a,nstep)call getv(v,a,nstep)enddodo 10 while(nstep.ne.0)print*,which step you want to know(100200),input 1 to endread*,nstepif(nstep.eq.0) stopcall displ(pos,v,nstep)10continueendc 设定初值subroutine begin(pos,v,a)dimension pos(64,3,200),v(64,3,200),a(64,3,200)do i=1,64do j=1,3pos(i,j,1)=int(rand()*10.0)v(i,j,1)=2.0*rand()-1.0 end doend doendc求得位置 velocity Verlet 算法subroutine getr(pos,v,a,nstep)dimension pos(64,3,200),v(64,3,200),a(64,3,200)parameter(dt=0.02)do i=1,64do j=1,3pos(i,j,nstep)=pos(i,j,nstep-1)+v(i,j,nstep-1)*dt $ -0.5D0*dt*dt*a(i,j,nstep-1) end doend doendc边界条件subroutine bound(pos,nstep)dimension pos(64,3,200),v(64,3,200),a(64,3,200)do i=1,64do j=1,3if(pos(i,j,nstep).le.0.0) thenpos(i,j,nstep)=pos(i,j,nstep)+10.0 elseif(pos(i,j,nstep).ge.10.0) thenpos(i,j,nstep)=pos(i,j,nstep)-10.0 end if end do end doendc 获得加速度subroutine getac(pos,a,nstep)dimension pos(64,3,200),v(64,3,200),a(64,3,200)double precision rij,xij,yij,zijforce(r)=48D0*r*(-13.0)-24D0*r*(-7.0)do i=1,64do j=i+1,64xij=pos(j,1,nstep)-pos(i,1,nstep)yij=pos(j,2,nstep)-pos(i,2,nstep)zij=pos(j,3,nstep)-pos(i,3,nstep)rij=sqrt(xij*2.0+yij*2.0+zij*2.0) if(rij.le.1E-5) rij=rij+0.02a(j,1,nstep)=a(j,1,nstep)+force(rij)/1.0*xij/rija(j,2,nstep)=a(j,2,nstep)+force(rij)/1.0*yij/rija(j,3,nstep)=a(j,3,nstep)+force(rij)/1.0*zij/rija(i,1,nstep)=a(i,1,nstep)-force(rij)/1.0*xij/rija(i,2,nst