应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

小木虫论坛-学术科研互动平台 » 计算模拟区 » 程序语言 » 【求助】请教fortran程序出错问题
51/1返回列表
查看: 514  |  回复: 4
只看楼主 @他人 存档 新回复提醒 (忽略) 收藏    在APP中查看
当前主题已经存档。

jianchaoyv

金虫 (小有名气)

[交流] 【求助】请教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
请多多指教！！
回复此楼

» 猜你喜欢
1楼 2009-04-13 16:31:04
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)
期待by 老虎大王的指点！！
赞一下 回复此楼
2楼2009-04-14 18:01:29
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
★ ★ ★ ★ ★
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()写在主程序的开头！！！呵呵。

注：我只管编程，没有考察算法的正确与否。
赞一下(1人) 回复此楼
3楼2009-04-14 19:45:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

jianchaoyv

金虫 (小有名气)
“第一个问题：为什么输出全是零呢？因为
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传到子程序里，当然输出的都是零。”
请问怎样把vx,vy,vz传到子程序里？？
赞一下 回复此楼
4楼2009-04-14 21:46:05
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

老虎大王

木虫 (著名写手)
★ ★
sunxiao(金币+2,VIP+0):谢谢参与，欢迎常来仿真编程版 4-15 09:00
啊？怎样传啊？比如你用个call comvel(temp, vx,vy,vz)吧。或者把vx,vy,vz  数组设为全局变量（用COMMON)。
赞一下(10人) 回复此楼
5楼2009-04-15 08:32:16
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 jianchaoyv 的主题更新
51/1返回列表
普通表情 龙 兔 虎 猫 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[找工作] 售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急 +3 YLHlRHNwYkce 2026-09-02 8/400 2026-09-03 12:05 by T0rGB46095mJ
[教师之家] 售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急 +3 YLHlRHNwYkce 2026-09-02 8/400 2026-09-03 12:03 by T0rGB46095mJ
[基金申请] 国社科又开始会评了，不知道这次命运如何 +9 雨打竹帘 2026-08-30 13/650 2026-09-03 08:16 by kudofaye
[基金申请] 麻烦专家们看看评委们的意见（F口面上） +8 gdd2018 2026-08-28 13/650 2026-09-03 06:22 by 羽毛枫f
[公派出国] 售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急 +3 ero8OE6tv9cu 2026-09-02 4/200 2026-09-03 03:14 by rM1TE0WVDIIY
[博后之家] 售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急 +3 ero8OE6tv9cu 2026-09-02 6/300 2026-09-02 23:59 by rM1TE0WVDIIY
[基金申请] 科研人应该花精力去思考如何解决问题，而不是去凝练问题 +7 瞬息宇宙 2026-09-01 14/700 2026-09-02 17:33 by ma0526
[基金申请] 要骂人了，新模版改版就是要淡化问题凝练这种虚的东西，结果有个评委还在说凝练得不够 +9 瞬息宇宙 2026-08-31 17/850 2026-09-02 14:04 by anjeeshine
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 10/500 2026-09-02 11:35 by 大不刘6
[基金申请] 学科评审组评审是指会评吗？ +5 瞬息宇宙 2026-08-31 5/250 2026-09-02 10:13 by 雨冰共舞
[基金申请] 面上函评意见出来了，像什么等级？ 20+4 Tsingking1 2026-08-27 17/850 2026-09-01 19:51 by 超级无敌华子
[论文投稿] 小白求助 投论文要求的highlights应该如何写 5+3 l1963982152 2026-08-29 4/200 2026-09-01 09:04 by 北京莱茵编辑
[基金申请] 面上意见出来了 +12 黄鸟于飞Chao 2026-08-29 23/1150 2026-08-31 18:57 by 黄鸟于飞Chao
[基金申请] 能否申诉？ +7 echo8914667 2026-08-30 8/400 2026-08-31 17:00 by yihongxu
[基金申请] 中青基了要发朋友圈吗？ +7 349506619 2026-08-28 7/350 2026-08-31 13:39 by 冼亮淀粉酶
[基金申请] 29号明天会评吗 +4 笨笨唐 2026-08-28 4/200 2026-08-31 09:30 by huixian257
[基金申请] 为什么到现在没收到通知？ +5 tannykie 2026-08-29 5/250 2026-08-30 21:05 by purplejack
[考博] 找导师 +6 yuanjiabao 2026-08-29 7/350 2026-08-30 14:40 by 生科新手
[基金申请] 我就是申请一个面上项目而已，这评审意见是按照杰青的条件评的吧？ +6 gouxfjh 2026-08-28 11/550 2026-08-30 07:57 by gouxfjh
[基金申请] 国自然评审意见 +13 wangmingqi 2026-08-28 19/950 2026-08-29 10:22 by Poppy1104
信息提示
关闭
请填处理意见
关闭