Fortran初學者 - 輸出文件顯示NaN

Question

我使用Euler-Richardson算法爲三體問題編寫了一個簡單的Fortran程序。出於某種原因，輸出文件只能提供NaN。這個問題是通過使用子程序還是函數來解決的？

PROGRAM Threebody 
IMPLICIT none 



!************** Variable declarations *************** 

!real(8), parameter :: G = 6.6738D-11 

INTEGER :: IOSTAT, dt, t_step, io_error, i 
DOUBLE PRECISION:: st, g, m0, m1, m2, force0x, force0y, force1x, force1y, force2x, force2y, n0, n1, n2, x0,v0,x1,v1,x2,v2,y0,w0,y1 
DOUBLE PRECISION:: w1,y2,w2 

m0 = 1.d0        ! the masses of the three bodies 
m1 = 1.d0 
m2 = 1.d0 

g = 9.80d0        ! m/s^2 s due to gravity 

t_step = 1 

y0 = 0.d0       
x0 = 0.d0 
v0 = 0.d0        
w0 = 0.d0 

y1 = 0.d0       
x1 = 1.d0 
v1 = 2.d0        
w1 = 1.118d0 

y2 = 0.d0       
x2 = -1.d0 
v2 = -1.118d0        
w2 = 0.d0 


dt = 1 ! time step 

OPEN(unit=5, file='out.txt', status='replace',action='write', IOSTAT=io_error) 
OPEN(unit=6, file='out1.txt', status='replace',action='write', IOSTAT=io_error) 
OPEN(unit=7, file='out2.txt', status='replace',action='write', IOSTAT=io_error) 





! particle no.1 
DO i=0,1000,t_step       ! 1000 is the total time 

st=i!/10.d0 

n0 = sqrt(((x2 - x1)**2 + (y2 - y1)**2)**3) 

force0x = (- m1*((x0 - x1)/n2)) - (m2*((x0 - x2)/n1)) 

x0 = x0 + (w0 + 0.5*force0x*st)*st 
w0 = w0 + force0x*st 

force0y = (- m1 * ((y0 - y1)/n2)) - (m2 * ((y0 - y2)/n1)) 

y0 = y0 + (v0 + 0.5*force0y*st)*st 
v0 = v0 + force0y*st 


WRITE(5,*) i*t_step, x0, y0, n0 

END DO 




! particle no.2 
DO i=0,1000,dt    

st=i 

n1 = sqrt(((x0 - x2)**2 + (y0 - y2)**2)**3) 

force1x = (- m2 * ((x1 - x0)/(n0))) - m0 * ((x1 - x2)/n2) 

x1 = x1 + (w1 + 0.5*force1x*st)*st 
w1 = w1 + force1x*st 

force1y = (- m2 * (y1 - y0)/n0) - (m0 * ((y1 - y2)/n2)) 

y1 = y1 + (v1 + 0.5*force1y*st)*st 
v1 = v1 + force1y*st 


WRITE(6,*) i*t_step, x1, y1 
END DO 




! particle no.3 
DO i=0,1000,dt 

st=i    

n2 = sqrt (((x1 - x0)**2 + (y1 - y0)**2)**3) 

force2x = (- m0 * ((x2 - x0)/n1)) - (m1 * ((x2 - x1)/n0)) 

x2 = x2 + (w2 + 0.5*force2x*st)*st 
w2 = w2 + force2x*st 

force2y = (- m0 * ((y2 - y0)/n1)) - (m1 * ((y2 - y1)/n0)) 

y2 = y2 + (v2 + 0.5*force2y*st)*st 
v2 = v2 + force2y*st 


WRITE(7,*) i*t_step, x2, y2, n2 
END DO 


CLOSE(unit=5) 
CLOSE(unit=6) 
CLOSE(unit=7) 

END PROGRAM Threebody 

Answer 1

在gfortran -Wall與編譯給你一個提示：

f.f90:56:0: warning: ‘n1’ may be used uninitialized in this function [-Wmaybe-uninitialized] 
     force0x = (- m1*((x0 - x1)/n2)) - (m2*((x0 - x2)/n1)) 

n1和n2不被初始化，使他們可以包含任何數據。在我的情況下，6.95e-310，導致強制給NaN。

要進行調試，您可以在屏幕上打印變量以檢查它們何時變爲NaN。爲此，請不要使用10以下的單位，因爲它們可能被保留。我認爲Unit 5是標準輸出，所以通過重新使用它你不能在屏幕上打印任何東西。