Fortran NaNの原因

Question

前提・実現したいこと

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)

! 初期条件
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)

! 彗星の位置推算　Runge-Kutta法

  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

! 微分方程式

  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

### 試したこと

sqrt内が正、分母に0が来ないことを確認しました。

### 補足情報（FW/ツールのバージョンなど）

数値計算で天体の軌道を求めるために運動方程式を解きたいと考えています。
Runge-Kutta法を用いるために運動方程式を1階の微分方程式に直して解く中で、
どうしても出力がNaNになってしまいます。また、その他の方程式は
f1=x4,f2=x5,f3=x6
となっております。
コンパイルはできており、またwrite文でバグの位置を探してみたところ、
関数(f4~f6)が怪しそうというところまでたどり着きましたが、心当たりがありません。
どなたか助言をいただけないでしょうか。

Answer 1

もし f4,f5,f6 に問題があるなら、たぶん

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

中の xj(i)*xj(i) が　xj(i)+xj(i) になっているところではと？　これが負数のため sqrt(負数) ⇒ NaN になっている可能性が高いと思います。。

これは最初に間違ったのをコピペしたようで３回繰り返して間違っています。

デバッグのヒントですが、NaN が出る場合を知っておくと少し楽になります。
0.0/0.0 （未定義演算）とか sqrt(負)　の様な組み込み関数の定義域未定義呼び出しとかが主なので、今の場合は割り算か sqrt が怪しいことになります。

Comments

[aaaaaaaaa_a]

ありがとうございます。完全に見落としていました。
NaNが出る場合の説明まで丁寧にありがとうございます。

Answer 2

xj,yj,zjの定義はどうなっているのでしょうか? 具体的な値が入っていないように見えますが…。

追記:

サブルーチン内のxy,yj,zjは、引数として渡すかcommon宣言を使うかしないかぎり、メインルーチンで定義されている同名の配列変数とは別物です。
いくらメイン側でreadしても、サブルーチン内の当該変数は未定義のままだと思います。

Comments

[aaaaaaaaa_a]

見てくださりありがとうございます。
定義は別のところでされているのですが、該当部分のプログラムだけを抜き出そうとしてかえって分かりにくくなってしまっていたのでソースコードを編集しました。
もしよければそちらを見ていただいて回答してくださると非常に助かります。

[aaaaaaaaa_a]

何度も対応してくださり本当にありがとうございました。

Answer 3

ありがとうございました。
(自己解決にしてしまっていました。申し訳ありません。答えてくださったお二方、ありがとうございました。)