Fortran Coder

查看: 185|回复: 5

[数值问题] fortran随机数调用问题

 关闭 [复制链接]

6

帖子

1

主题

0

精华

入门

F 币
27 元
贡献
25 点

规矩勋章

发表于 2017-5-11 08:41:09 | 显示全部楼层 |阅读模式
20F 币
本帖最后由 老衲不吸烟 于 2017-5-11 09:41 编辑

c(i)=normal(5.0d0,2.0d0)为什么没用啊,请教大神指条明路

回复

使用道具 举报

6

帖子

1

主题

0

精华

入门

F 币
27 元
贡献
25 点

规矩勋章

 楼主| 发表于 2017-5-11 08:42:50 | 显示全部楼层
[Fortran] 纯文本查看 复制代码
module phase
  use imslf90
  implicit none
  real*8,parameter,private:: pi=3.141592653589793238462643d0
  integer,parameter,private:: N=233, NN=N*2, M=1, tran=2000, dot=2000, npara=101,D=0.001d0
  !integer,parameter,private:: mxparm=50, N=30, NN=N*2, M=1, tran=200, dot=200, npara=101 !Fig. 4
  real*8, parameter,private:: gamma=0.7d0, a=1.d0, b=144.d0/233.d0, beta=89.d0/144.d0, L=144.d0, Kappa=3.d0*pi
  !real*8, parameter,private:: gamma=0.7d0, a=1.d0, b=1.d0, beta=24.d0/30.d0, L=30.d0, Kappa=3.d0*pi !Fig. 4
  real*8, private:: changeF
  private derivs
  public plot 
  contains

  function ran()   !returns random number between 0 - 1  
    implicit none 
    integer , save :: flag = 0
    double precision :: ran 
    if(flag==0) then 
      call random_seed()
      flag = 1 
    endif 
    call random_number(ran)     ! built in fortran 90 random number function  
  end function ran
  
  function normal(mean,sigma) 
    implicit none 
    integer :: flag 
    double precision, parameter :: pi = 3.141592653589793239  
    double precision :: uuu1, uuu2, yyy1, yyy2, normal, mean, sigma 
    save flag 
    data flag /0/ 
    uuu1 = ran(); uuu2 = ran() 
    if (flag.eq.0) then 
      yyy1 = sqrt(-2.0d0*log(uuu1))*cos(2.0d0*pi*uuu2) 
      normal = mean + sigma*yyy1 
      flag = 1 
    else 
      yyy2 = sqrt(-2.0d0*log(uuu1))*sin(2.0d0*pi*uuu2) 
      normal = mean + sigma*yyy2 
      flag = 0 
    endif  
  end function normal
 
 subroutine plot()
  integer i,j,k
  real*8 y(NN),dydt(NN)
  real*8 hh,th,h,t
  real*8 dym(NN),dyt(NN),yt(NN),c(N)
  real*8 y0(NN)
  real*8 parameters(npara),averageVelocity(npara),parameters2(npara),averageVelocity2(npara),temp
  character*38 filenm
  ! set initial conditions
  t = 0.0d0
  !call random_seed()
  do i=1,N
    !call random_number(y(i))
    y0(i)  =(i-1)*b
	y0(i+N)=0.d0
  enddo


!绝热增加F
  do k=1,npara
  y=y0
  changeF= 0.01d0*(k-1)
  t=0.d0; h=0.001d0
  hh=h/2.d0
  temp=0.d0
  do i=1,dot
     call derivs(t,y,dydt)
     yt=y+h*dydt
	 t=t+h
     call derivs(t,yt,dyt) 
     y=y+hh*(dydt+dyt) 
     do j=1,N
        temp=temp+y(j+N)
     enddo
  enddo
  temp=temp/N/dot
  print*, k,changeF,temp
  averageVelocity(k)=temp
  parameters(k)=changeF
  y0=y
  enddo


!输出结果
  filenm="eeee.txt"
  open (5555, file = filenm)
  write (5555,*) " x= [ "
  do i = 1 , npara 
     write(5555,'(f13.5)') parameters(i)
  enddo
  write (5555,*) " ]; "
  write (5555,*) " y= [ "
  do i = 1 , npara 
     write(5555,'(f13.5)') averageVelocity(i) 
  enddo
  write (5555,*) " ]; "
  write (5555,*) 
  write (5555,*) "figure(8)"
  write (5555,*) "plot(x,y,'-b^');"
  print*,filenm
  close (5555)
 endsubroutine plot

subroutine derivs(t,y,dydt)
implicit none
real*8 t
real*8 y(NN),dydt(NN),c(N)
integer i,j
real*8 u(N),v(N),w(-N+1:NN),PForce,f1,f2
  do i=1,N
     u(i)=y(i); v(i)=y(i+N)
	 c(i)=normal(5.0d0,2.0d0)
	 !为了处理周期性边界条件Eq.(3),引入如下中间变量
	 w(i)=y(i); w(i+N)=y(i)+N*b; w(i-N)=y(i)-N*b
  enddo
  do i=1,N
     PForce=0.d0;f1=0.d0;f2=0.d0
     do j=1,M
	    f1=f1-Kappa*(1.d0-dexp(b+w(i)-w(i+j))) * dexp(b+w(i)-w(i+j))   ! when x_i<x_j
		f2=f2+Kappa*(1.d0-dexp(b-w(i)+w(i-j))) * dexp(b-w(i)+w(i-j))
        PForce=f1+f2
	 enddo
     dydt(i)  =v(i)+D*c(i)
	 dydt(i+N)=-gamma*v(i) - dsin(2.d0*pi/a*u(i))/2.d0 - dsin(2.d0*pi*beta/a*u(i))/2.d0 -PForce+changeF+D*c(i)
  enddo
end subroutine  derivs

endmodule phase



program phase_plot
  use phase
  implicit none
  real*8 t1,t2
  call cpu_time(t1)
  call plot()
  call cpu_time(t2)
  print*, "Time used:", t2-t1, "seconds."
endprogram phase_plot

回复

使用道具 举报

6

帖子

1

主题

0

精华

入门

F 币
27 元
贡献
25 点

规矩勋章

 楼主| 发表于 2017-5-11 08:43:50 | 显示全部楼层
完整代码在一楼,求大神指出问题。。。
回复

使用道具 举报

6

帖子

1

主题

0

精华

入门

F 币
27 元
贡献
25 点

规矩勋章

 楼主| 发表于 2017-5-11 09:40:27 | 显示全部楼层
不能删帖么,逗比了,D=0了。。。
回复

使用道具 举报

6

帖子

1

主题

0

精华

入门

F 币
27 元
贡献
25 点

规矩勋章

 楼主| 发表于 2017-5-11 09:41:46 | 显示全部楼层
本帖最后由 老衲不吸烟 于 2017-5-11 09:43 编辑

ok。。。。。。。。。。。
回复

使用道具 举报

1096

帖子

12

主题

5

精华

论坛跑堂

Fcode跑堂

F 币
847 元
贡献
746 点

新人勋章贡献勋章管理勋章帅哥勋章爱心勋章规矩勋章元老勋章水王勋章

发表于 2017-5-11 18:52:22 | 显示全部楼层
没有开放删帖的功能。我可以帮你关闭帖子。
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|QQ群|Fcode

GMT+8, 2017-8-22 03:35

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表