計(jì)算物理實(shí)驗(yàn)考試及答案.doc

07級(jí)計(jì)算物理實(shí)驗(yàn)期末考試試題解答姓名：馬瑞松 學(xué)號(hào)：020112007016 班級(jí)：物理學(xué)2007級(jí) 成績(jī)： 1、 用fortran 語言編程采用兩種方法（高等數(shù)學(xué)方法和montecarlo方法）計(jì)算。答： 1、高等數(shù)學(xué)方法：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 此程序經(jīng)調(diào)試無誤，輸出結(jié)果為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此程序輸出結(jié)果為0.68860，以mtkljifen.f保存在所提交所提交壓縮包里。2、 用Gear predictor-corrector算法數(shù)值求解正方形區(qū)域 的拉普拉斯方程答：所編程序如下所示：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此程序經(jīng)過調(diào)試后以laplace.f保存在所提交所提交壓縮包里。3、 十個(gè)氫原子在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個(gè)符合麥克斯韋分布的速度，此程序以maikesiwei.f保存在所提交壓縮包內(nèi)。4、 十個(gè)氫原子在模擬盒中運(yùn)動(dòng)，如果盒中有1000米/秒得風(fēng)自左向右吹，用fortran語言編程計(jì)算左右器壁的壓力差，適當(dāng)選取模擬盒大小和溫度，原子與器壁的碰撞是完全彈性的。答：所編程序如下所示： 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的正方體，由鍵盤輸入氫原子相對(duì)原子質(zhì)量1以及自行設(shè)定的溫度，其初速度符合麥克斯韋速率分布，輸出結(jié)果為左右器壁的壓強(qiáng)差與壓力差。調(diào)試后程序以press.f保存在所提交壓縮包內(nèi)。五、編寫一個(gè)三維，元胞尺寸為L3的周期邊界條件計(jì)算程序。答：所編程序如下所示：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此程序假設(shè)元胞長度L為1，模擬一個(gè)原子在元胞內(nèi)采用周期性邊界條件的運(yùn)動(dòng)。原子的初始位置以及初始速度隨機(jī)給出，由鍵盤輸入欲想知道的步數(shù)，即可得到該原子此時(shí)的位置以及周期性邊界條件采用的次數(shù)。此程序以boundary.f保存在所提交壓縮包內(nèi)。6、 試做總能量固定的單原子系統(tǒng)的分子動(dòng)力學(xué)模擬。元胞為，劃分為的正方形網(wǎng)格。元胞內(nèi)原子數(shù)。原子質(zhì)量。位勢(shì)為Lenard-Jones勢(shì)，其中，邊界條件為周期性邊界條件，初始位置是隨機(jī)分布在正則節(jié)點(diǎn)上，初始速度為按-1，1隨機(jī)分布。分子動(dòng)力學(xué)模擬步長取為，模擬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 設(shè)定初值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