fortran 代码不循环

Question

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

           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，有人可以帮帮我吗？

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

Answer 1

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

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

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

其他地方也一样