Home>

I want to remove the cause of NaN.

```
###
program main
implicit none
integer,parameter :: n=14610
integer :: i,day,k
real*8 :: h,a(n),xj(n),yj(n),zj(n),x1,x2,x3,x4,x5,x6,
&f1,f2,f3,f4,f5,f6,ka1,ka2,ka3,ka4,ka5,ka6,kb1,kb2,kb3,
&kb4,kb5,kb6,kc1,kc2,kc3,kc4,kc5,kc6,kd1,kd2,kd3,kd4,
&kd5,kd6,gs,gj,xr,xrj,xmj
open(10,file='z_jupitar-2.txt')
do i=1,n
read(10,*) a(i),xj(i),yj(i),zj(i)
end do
close(10)
!Initial conditions
x1=3.64d-7
x2=0.0d0
x3=0.0d0
x4=0.0d0
x5=8.0d0
x6=0.0d0
h=20.0d0
write(*,*) 0,x1,x2,x3,xj(1),yj(1),zj(1)
!Comet position estimation Runge-Kutta method
do k=1,n-3,2
i=k
day=(i+1)*10
ka1=h*f1(x4)
ka2=h*f2(x5)
ka3=h*f3(x6)
ka4=h*f4(x1,x2,x3)
ka5=h*f5(x1,x2,x3)
ka6=h*f6(x1,x2,x3)
!write(*,*) ka1,ka2,ka3,ka4,ka5,ka6
i=i+1
kb1=h*f1(x4+ka4/2.0d0)
kb2=h*f2(x5+ka5/2.0d0)
kb3=h*f3(x6+ka6/2.0d0)
kb4=h*f4(x1+ka1/2.0d0,x2+ka2/2.0d0,x3+ka3/2.0d0)
kb5=h*f5(x1+ka1/2.0d0,x2+ka2/2.0d0,x3+ka3/2.0d0)
kb6=h*f6(x1+ka1/2.0d0,x2+ka2/2.0d0,x3+ka3/2.0d0)
!write(*,*) kb1,kb2,kb3,kb4,kb5,kb6
kc1=h*f1(x4+kb4/2.0d0)
kc2=h*f2(x5+kb5/2.0d0)
kc3=h*f3(x6+kb6/2.0d0)
kc4=h*f4(x1+kb1/2.0d0,x2+kb2/2.0d0,x3+kb3/2.0d0)
kc5=h*f5(x1+kb1/2.0d0,x2+kb2/2.0d0,x3+kb3/2.0d0)
kc6=h*f6(x1+kb1/2.0d0,x2+kb2/2.0d0,x3+kb3/2.0d0)
!write(*,*) kc1,kc2,kc3,kc4,kc5,kc6
i=i+1
!write(*,*) i
kd1=h*f1(x4+kc4)
kd2=h*f2(x5+kc5)
kd3=h*f3(x6+kc6)
kd4=h*f4(x1+kc1,x2+kc2,x3+kc3)
kd5=h*f5(x1+kc1,x2+kc2,x3+kc3)
kd6=h*f6(x1+kc1,x2+kc2,x3+kc3)
!write(*,*) kd1,kd2,kd3,kd4,kd5,kd6
x1=x1+(ka1+kd1)/6.0d0+(kb1+kc1)/3.0d0
x2=x2+(ka2+kd2)/6.0d0+(kb2+kc2)/3.0d0
x3=x3+(ka3+kd3)/6.0d0+(kb3+kc3)/3.0d0
x4=x4+(ka4+kd4)/6.0d0+(kb4+kc4)/3.0d0
x5=x5+(ka5+kd5)/6.0d0+(kb5+kc5)/3.0d0
x6=x6+(ka6+kd6)/6.0d0+(kb6+kc6)/3.0d0
write(*,*) day,x1,x2,x3,x4,x5,x6
end do
stop
end program main
Differential equation
real*8 function f1(x4)
implicit none
real*8 :: x4
f1=x4
return
end function f1
real*8 function f2(x5)
implicit none
real*8 :: x5
f2=x5
return
end function f2
real*8 function f3(x6)
implicit none
real*8 :: x6
f3=x6
return
end function f3
real*8 function f4(x1,x2,x3)
implicit none
integer,parameter :: n=14610
integer :: i
real*8 :: x1,x2,x3,xj(n),yj(n),zj(n),
&gs,gj,xr,xrj,xmj
gs=2.95912208d-4
gj=2.82534589d-7
xr=sqrt(x1*x1+x2*x2+x3*x3)
xrj=sqrt(xj(i)+xj(i)+yj(i)*yj(i)+zj(i)*zj(i))
xmj=sqrt((xj(i)-x1)*(xj(i)-x1)+(yj(i)-x2)*(yj(i)-x2)
&+(zj(i)-x3)*(zj(i)-x3))
f4=-gs*x1/(xr*xr*xr)
&+gj*((xj(i)-x1)/(xmj*xmj*xmj)-xj(i)/(xrj*xrj*xrj))
return
end function f4
real*8 function f5(x1,x2,x3)
implicit none
integer,parameter :: n=14610
integer :: i
real*8 :: x1,x2,x3,xj(n),yj(n),zj(n),
&gs,gj,xr,xrj,xmj
gs=2.95912208d-4
gj=2.82534589d-7
xr=sqrt(x1*x1+x2*x2+x3*x3)
xrj=sqrt(xj(i)+xj(i)+yj(i)*yj(i)+zj(i)*zj(i))
xmj=sqrt((xj(i)-x1)*(xj(i)-x1)+(yj(i)-x2)*(yj(i)-x2)
&+(zj(i)-x3)*(zj(i)-x3))
f5=-gs*x2/(xr*xr*xr)
&+gj*((yj(i)-x2)/(xmj*xmj*xmj)-yj(i)/(xrj*xrj*xrj))
return
end function f5
real*8 function f6(x1,x2,x3)
implicit none
integer,parameter :: n=14610
integer :: i
real*8 :: x1,x2,x3,xj(n),yj(n),zj(n),
&gs,gj,xr,xrj,xmj
gs=2.95912208d-4
gj=2.82534589d-7
xr=sqrt(x1*x1+x2*x2+x3*x3)
xrj=sqrt(xj(i)+xj(i)+yj(i)*yj(i)+zj(i)*zj(i))
xmj=sqrt((xj(i)-x1)*(xj(i)-x1)+(yj(i)-x2)*(yj(i)-x2)
&+(zj(i)-x3)*(zj(i)-x3))
f6=-gs*x3/(xr*xr*xr)
&+gs*((zj(i)-x3)/(xmj*xmj*xmj)-zj(i)/(xrj*xrj*xrj))
return
end function f6
What I tried

I confirmed that sqrt is positive and that 0 does not come in the denominator.

Supplementary information (FW/tool ​​version, etc.)

I would like to solve the equation of motion to obtain the orbit of the celestial body by numerical calculation.
In order to use the Runge-Kutta method by solving the equation of motion by solving it into a first-order differential equation,
The output will always be NaN. And other equations are
f1=x4,f2=x5,f3=x6
It has become.
It is compiled, and when I searched for the position of the bug with the write statement,
I reached the point where the functions (f4 to f6) seemed suspicious, but I have no idea.
Would anyone please give me some advice?

  • Answer # 1

    What's the definition of xj, yj, zj? It doesn't seem to contain any specific values...

    Addendum:

    The xy,yj,zj in the subroutine are different from the array variables of the same name defined in the main routine unless you pass them as arguments or use the common declaration.
    No matter how much you read on the main side, I think that the variable in the subroutine remains undefined.

  • Answer # 2

    If you have problems with f4,f5,f6, maybe

    xrj=sqrt(xj(i)+xj(i)+yj(i)*yj(i)+zj(i)*zj(i))


    Where is xj(i)*xj(i) inside is xj(i)+xj(i)? Since this is a negative number, it is highly possible that sqrt (negative number) ⇒ NaN. ..

    It seems that I copied and copied the wrong one at the beginning, so I repeated it three times and made a mistake.

    It's a debugging tip, but it's a little easier to know when you get NaN.
    Since it is mainly the domain undefined call of built-in functions such as 0.0/0.0 (undefined operation) or sqrt (negative), division or sqrt is questionable in this case.

  • Answer # 3

    https://stackoverflow.com/questions/276964?modal=q-comp

    I have a question following.
    I've rewritten the code a lot, so I'll solve this question for now.
    Thank you all.