老衲不吸烟 发表于 2017-5-11 08:41:09

fortran随机数调用问题

本帖最后由 老衲不吸烟 于 2017-5-11 09:41 编辑

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

老衲不吸烟 发表于 2017-5-11 08:42:50

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 subroutinederivs

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

老衲不吸烟 发表于 2017-5-11 08:43:50

完整代码在一楼，求大神指出问题。。。

老衲不吸烟 发表于 2017-5-11 09:40:27

不能删帖么，逗比了，D=0了。。。

老衲不吸烟 发表于 2017-5-11 09:41:46

本帖最后由 老衲不吸烟 于 2017-5-11 09:43 编辑

ok。。。。。。。。。。。

fcode 发表于 2017-5-11 18:52:22

没有开放删帖的功能。我可以帮你关闭帖子。
页: [1]
查看完整版本: fortran随机数调用问题