| 查看: 428 | 回复: 4 | |||
| 当前主题已经存档。 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[交流]
【求助】请教fortran程序出错问题
|
|||
|
问题一：程序如下： program initial_velocity implicit none !we assumed m=1,k=1 real*8,parameter::temp=300.0 integer,parameter::N=10 integer::i real*8,dimension(1:N)::vx,vy,vz write(*,'(1x,a,8x,a,10x,a,10x,a)')'n','vx','vy','vz' call comvel(temp) do i=1,N write(*,'(i2,2x,f10.3,2x,f10.3,2x,f10.3)')i,vx(i),vy(i),vz(i) end do end program initial_velocity subroutine comvel(temp) integer,parameter::N=10 integer::i real*8,dimension(1:N)::vx,vy,vz real*8::temp,Rtemp,sumx,sumy,sumz real*8::gauss,dummy Rtemp=sqrt(temp) do i=1,N vx(i)=Rtemp*gauss(dummy) vy(i)=Rtemp*gauss(dummy) vz(i)=Rtemp*gauss(dummy) end do sumx=0.0 sumy=0.0 sumz=0.0 do i=1,N sumx=sumx+vx(i) sumy=sumy+vy(i) sumz=sumz+vz(i) end do sumx=sumx/real(N) sumy=sumy/real(N) sumz=sumz/real(N) do i=1,N vx(i)=vx(i)-sumx vy(i)=vy(i)-sumy vz(i)=vz(i)-sumz end do return end subroutine comvel function gauss(dummy) real*8,parameter::A1=3.949846138,A3=0.252408784 real*8,parameter::A5=0.076542912,A7=0.008355968 real*8,parameter::A9=0.029899776 real*8::sum,R,R2,dummy integer::i sum=0.0 call random_seed() do i=1,12 call random_number(dummy) sum=sum+dummy end do R=(sum-6.0)/4.0 R2=R*R gauss=((((A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R return end function gauss 运行后结果： n vx vy vz 1 0.000 0.000 0.000 2 0.000 0.000 0.000 3 0.000 0.000 0.000 4 0.000 0.000 0.000 5 0.000 0.000 0.000 6 0.000 0.000 0.000 7 0.000 0.000 0.000 8 0.000 0.000 0.000 9 0.000 0.000 0.000 10 0.000 0.000 0.000 Press any key to continue 请问怎样能得到所需的随机数值？？ 问题二：为验证上面程序中部分程序段，我写另一程序如下： program a_1 implicit none real*8::gauss,dummy,Rtemp integer,parameter::N=10 integer::i real,dimension(1:N)::vx,vy,vz real,parameter::temp=300.0 Rtemp=sqrt(temp) do i=1,N vx(i)=Rtemp*gauss(dummy) vy(i)=Rtemp*gauss(dummy) vz(i)=Rtemp*gauss(dummy) print*,vx(i),vy(i),vz(i) end do end program a_1 function gauss(dummy) real*8,parameter::A1=3.949846138,A3=0.252408784 real*8,parameter::A5=0.076542912,A7=0.008355968 real*8,parameter::A9=0.029899776 real*8::sum,R,R2,dummy integer::i sum=0.0 do i=1,12 call random_seed() call random_number(dummy) sum=sum+dummy end do R=(sum-6.0)/4.0 R2=R*R gauss=((((A9*R2+A7)*R2+A5)*R2+A3)*R2+A1)*R return end function gauss 运行结果有时值全部一样，有时部分一样，但对应同一个i 时vx(i),vy(i),vz(i)始终值一样结果如下： 结果（1）： 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 49.20191 Press any key to continue 结果（2）： 82.36326 82.36326 82.36326 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 -72.15305 Press any key to continue 请多多指教！！ |
回复此楼» 猜你喜欢
一志愿华东师范大学有机化学专业，初试351分，复试被刷求调剂!
已经有4人回复
-
298求调剂
已经有4人回复
-
材料学硕333求调剂
已经有11人回复
-
0856求调剂
已经有7人回复
-
356求调剂
已经有3人回复
-
332求调剂
已经有7人回复
-
348求调剂
已经有4人回复
-
305求调剂
已经有8人回复
-
一志愿北京化工大学材料与化工（085600）296求调剂
已经有13人回复
-
调剂310
已经有9人回复
4楼2009-04-14 21:46:05
2楼2009-04-14 18:01:29
老虎大王
木虫 (著名写手)
- 应助: 26 (小学生)
- 贵宾: 0.17
- 金币: 4774.1
- 散金: 8
- 红花: 42
- 帖子: 1361
- 在线: 215.2小时
- 虫号: 659094
- 注册: 2008-11-21
- 专业: 金属结构材料
★ ★ ★ ★ ★
jianchaoyv(金币+5,VIP+0): 4-14 21:25
jianchaoyv(金币+5,VIP+0): 4-14 21:25
|
第一个问题：为什么输出全是零呢？因为 call comvel(temp) do i=1,N write(*,'(i2,2x,f10.3,2x,f10.3,2x,f10.3)')i,vx(i),vy(i),vz(i) end do 在这里，你没有把vx,vy,vz传到子程序里，当然输出的都是零。 第二个问题，我建议你运行以下这两个小程序，测试一下 random_seed()和random_number(dummy)的用法。 程序一：（这是你的程序写法） do i=1,12 call random_seed() call random_number(dummy) write(*,*) dummy enddo end 得到什么？呵呵，全都一样的随机数。但是，每运行一次，就有不同的序列。 程序二： 把random_seed()放在外面： call random_seed() do i=1,12 call random_number(dummy) write(*,*) dummy enddo end 得到什么？不同的随机数！ 所以，你应该把 call random_seed()写在循环的外边！ 但是，你把Call random_seed()写在子程序中循环的外边之后，又如何呢？ 又是同样的随机数！！！为什么呢？因为你在主程序中反复调用了Gauss，还是相当于把Call random_seed()写在了循环内部。 最后我们知道，正确的写法应该把Call random_seed()写在主程序的开头！！！呵呵。 注：我只管编程，没有考察算法的正确与否。 |
3楼2009-04-14 19:45:56
老虎大王
木虫 (著名写手)
- 应助: 26 (小学生)
- 贵宾: 0.17
- 金币: 4774.1
- 散金: 8
- 红花: 42
- 帖子: 1361
- 在线: 215.2小时
- 虫号: 659094
- 注册: 2008-11-21
- 专业: 金属结构材料
5楼2009-04-15 08:32:16
