fortran 代码不循环

Posted 2023-03-11

【中文标题】fortran 代码不循环

【英文标题】：fortran code do not cycle

【发布时间】：2014-02-25 14:39:44

【问题描述】：

我尝试用以下代码解决这个简单的三体问题：

           Program Main
           Implicit real*8 (A-H,O-Z)
           real*8 ome,mu, rho, R
           duepi=8*datan(1.d0)
           ome=1
           mu=0.001


           T_per=duepi/ome


           rho=0.1
           R=1.0
           N_step=100

c Open the file              
             OPEN(unit=11, file="prova1.txt")
c Nested do loops
             do iy0=1,100
               do iP0=1,100
c Calc value for 0
                 y0 = real(iy0)/100.
                 x0 = 0
c Calc value for py0
                 py0 = real(iP0)/100.
                 px0 = 0

           x=x0
           y=y0
           px=px0
           py=py0


           dt=T_per/N_step

           E0=H(x,y,px,py)

           k_max=100*N_step

           k=0
           t=0
           errh=0

c---------
c  start integration loop
c--------
           do k=1,k_max
           call sym4(x,y,px,py,dt)
           E= H(x,y,px,py)
           errh=abs(E-E0)
           t=k*dt

           enddo



           do k=1,k_max
           call sym4(x,y,px,py,-dt)
           E= H(x,y,px,py)
           errh=abs(E-E0)
           t=t-dt

           enddo

        write(11,*) y0, py0, errh 
               enddo ! iP0
             enddo ! iy0
        close(11)

           end



           subroutine sym1(x,y,px,py,dt)
           Implicit real*8 (A-H,O-Z)
c
           call f(x,y,fx,fy)

           pxnew=px+dt*fx
           pynew=py+dt*fy


           xnew=x+dt*pxnew
           ynew=y+dt
c
           x=xnew
           y=ynew
           px=pxnew
           py=pynew

           end

           subroutine sym1_B(x,y,px,py,dt)
           Implicit real*8 (A-H,O-Z)
c

           xnew=x+dt*px
           ynew=y+dt

           call f(xnew,ynew,fxnew,fynew)
           pxnew=px+dt*fxnew
           pynew=py+dt*fynew 
c
           x=xnew
           y=ynew
           px=pxnew
           py=pynew

           end



           subroutine  f(x,y,fx,fy)
           Implicit real*8 (A-H,O-Z)
           real*8 ome,mu,rho,R


            fx = ((1-mu)*(rho+x))/((rho*rho+2*rho*x+y*y)**(1.5)) -
     &           (mu*(R+x))/((R**2-2*R*x+x*x+y*y)**(1.5))
            fy = ((1-mu)*(rho+x))/((rho**2+x**2+2*rho*x+y**2)**(1.5))/
     &       + (mu*y)/((R**2-2*R*x+x**2+y**2)**(1.5))   

           return
           end

           real*8 function H(x,y,px,py)
           Implicit real*8 (A-H,O-Z)
           real*8 ome,mu,rho,R

c          h=px*px/2.d0+ py +(1+eps*cos(ome*y))*x*x/2

c           h=px*px/2.d0+ py -(1+eps*cos(ome*y))*cos(x) 

    r12 = sqrt( ( (x*cos(ome*y)-y*sin(ome*y))+rho*cos(ome*y) )**2
     &        + ( x*sin(ome*y)+y*cos(ome*y) + rho*sin(ome*y) )**2 )
        r13 = sqrt( ( (x*cos(ome*y)-y*sin(ome*y))-R*cos(ome*y) )**2
     &        + ( x*sin(ome*y)+y*cos(ome*y) - R*sin(ome*y) )**2 )       

        h=(px**2+py**2)/2.d0 - (1-mu)/r12 - mu/r13 + py

           return
           end


           subroutine sym2(x,y,px,py,dt)
           Implicit real*8 (A-H,O-Z)

           call f(x,y,fx,fy)


           xnew= x+ px*dt +    fx*dt**2/2.d0   
           ynew= y+ dt                         ! così   è giusto

           call f(xnew,ynew,fxnew,fynew) 
           pxnew= px+ dt*(fx+fxnew )/2.d0
           pynew= py+ dt*(fy+fynew )/2.d0

           x=xnew
           y=ynew
           px=pxnew
           py=pynew

           end


           subroutine sym4(x,y,px,py,dt)
           Implicit real*8 (A-H,O-Z)
           sq2=2**(1.d0/3.d0)
           alpha= 1.d0/(2-sq2)
           beta= sq2/(2-sq2)
           dt1= dt*alpha
           dt2=-dt*beta
           call sym2(x,y,px,py,dt1)
           call sym2(x,y,px,py,dt2)
           call sym2(x,y,px,py,dt1)
           return
       end

代码调用辛积分器并解决 3body 问题。但是当我尝试运行它时，没有编译错误，output.txt 文件只显示初始网格而不是errh 这个专栏只给我 NaN，有人可以帮帮我吗？

可能是初始条件错误（速度或位置的奇怪初始条件，欧米茄是错误的......）？

【问题讨论】：

这里不足以理解问题。我建议添加一些中间写语句。最重要的是，我建议使用 Fortran 95 而不是 FORTRAN 77。

检查你为什么得到错误结果的常用方法是检查算法，也许使用调试器。

任何在 C21 中编写 Fortran 而不在所有范围内使用 implicit none 的人都是自找麻烦。隐式类型是否是您当前问题的根源并不重要。重要的是你，任何人，都应该使用你的语言和编译器提供的所有工具来消除错误的可能性。就个人而言，我什至不会尝试在没有implicit none 的情况下调试代码。

【参考方案1】：

正如 cmets 中所述，您的程序存在许多代码样式问题（意图错误，...）以及未使用最佳 Fortran 实践（例如 implicit none）。您的问题可以通过简单地使用调试器来解决。

一个明显的问题是您在函数中使用了未初始化的变量：

subroutine f(x,y,fx,fy) : rho, mu, R 在计算 fx 和 fy，产生 NaN function H()：ome, rho, mu，...

其他地方也一样