Fortran Coder

标题: 负数求三次方根的NaN问题，多谢啦！ [打印本页]

作者: 小智 时间: 2015-9-30 01:42
标题: 负数求三次方根的NaN问题，多谢啦！
不说多的直接上源代码

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode



program ex110

implicit none

DIMENSION U(2),UI1(3),UI2(6),UI3(6) 

REAL(8) Nv, chi, lambda0, detF0,vcho,vcmo,vcpo,fa,fb,Ka,Kb,vNa

REAL(8) U ,UI1, UI2 ,UI3,aa, bb, cc ,dd, ee, deltaa, deltab,deltac

REAL(8) BI1,BI2,AJ

REAL(8) delta,vch,GAJ, DGAJDAJ, DDGAJDDAJ, DDDGAJDDDAJ 

REAL(8) DGAJDvch ,DDGAJDDvch, DDDGAJDDDvch 

REAL(8) DDGAJDAJDvch, DDDGAJDAJDDvch, DDDGAJDDAJDvch

REAL(8) DvchDAJ, DDvchDDAJ, DDDvchDDDAJ

REAL(8) D1GAJDAJ, D2GAJDAJ, D3GAJDAJ,TEMP

BI1=10.0

AJ=1.0

TEMP=2.0

Nv = 0.001

chi = 0.1

lambda0 = 3.39

detF0 = lambda0**3.0

vNa=0.0602

vcpo=vNa*(10.0**(-3.0))

vcho=vNa*(10.0**(-TEMP))

vcmo=vcpo+vcho

Ka=10.0**(-4.3)

Kb=10.0**(-10.0)

fa=0.05

fb=0.0

!。。。。。。中间省略部分无关代码





aa=AJ*detF0*(vcho+vcpo)



bb=fb*vcho+AJ*detF0*Ka*(vcho+vcpo)*vNa+ AJ*detF0*Kb*(vcho + 

&vcpo)*vNa



cc=-(AJ*detF0*vcho**2*vcmo)-fa*Ka*vcho*vNa+fb*Ka*vcho*vNa+ 

&AJ*detF0*Ka*Kb*(vcho+vcpo)*vNa**2

dd=-(AJ*detF0*Ka*vcho**2*vcmo*vNa)-AJ*detF0*Kb*vcho**2*vcmo*v

&Na-fa*Ka*Kb*vcho*vNa**2



ee=AJ*detF0*Ka*Kb*vcho**2*vcmo*vNa**2



deltaa=cc**2-3*bb*dd+12*aa*ee



deltab=2*cc**3-9*bb*cc*dd+27*aa*dd**2-72*aa*cc*ee

deltac=(deltab+Sqrt(abs(-4.0*deltaa**3.0+deltab**2.0)))

&**(1.0/3.0)



delta=(2**(1/3)*deltaa)/(3*aa*deltac)+deltac/(3*2**(1/3)*aa)

!。。。。。。。。。。。以下省略部分无关代码。。。。。。。。。

现在我的问题是，deltac=(deltab+Sqrt(abs(-4.0*deltaa**3.0+deltab**2.0)))
&**(1.0/3.0)这一句中，如果只有(deltab+Sqrt(abs(-4.0*deltaa**3.0+deltab**2.0)))的话能运行出结果来
但是如果我想将这个deltac开立方根，变成deltac=(deltab+Sqrt(abs(-4.0*deltaa**3.0+deltab**2.0)))**(1.0/3.0)之后
运行结果死活就是NaN，实在不知道是为什么求大神指点多谢~

PS：打印出来的deltac就是此图中第七行第二个数据，如果没有开立方根的话能运行出结果，开了立方根之后死活是NaN。。。。。T_T

作者: 楚香饭 时间: 2015-9-30 07:35
本帖最后由楚香饭于 2015-9-30 07:36 编辑

由于浮点数有误差。所以 1.0/3.0 并不能精确的表示 1/3，例如可能是 0.333333334，极大可能是一个偶数分之一。
你知道，负数开偶次方会是 NaN 对吧？

所以，如果要开三次方，比较好的办法，是先取绝对值，开方以后，再把符号乘上去。（sign函数可以返回符号，即1.0或-1.0）

以下代码在三个编译器上测试，均能得到正确结果。

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

Program www_fcode_cn

  Implicit None

  real(8) :: d = -1.63899782385934d-24 , s

  s = sign(1.0d0,d)*abs(d)**(1.0d0/3.0d0)

  write(*,*) s , d , s*s*s

End Program www_fcode_cn

作者: 小智 时间: 2015-9-30 08:30

楚香饭发表于 2015-9-30 07:35
由于浮点数有误差。所以 1.0/3.0 并不能精确的表示 1/3，例如可能是 0.333333334，极大可能是一个偶数分之 ...

果然编译成功了谢谢大神！

作者: 糖盒love玲珑 时间: 2015-12-28 19:18
楚兄这个例子好可以作为新手认识浮点数误差的一个很好案例! 受教了.

作者: pasuka 时间: 2016-1-2 20:14

楚香饭发表于 2015-9-30 07:35
由于浮点数有误差。所以 1.0/3.0 并不能精确的表示 1/3，例如可能是 0.333333334，极大可能是一个偶数分之 ...

查了一下，C++11的标准里面已经有一个cbrt函数专门开三次方根，不知道Fortran何时跟进

作者: wutong-sky 时间: 2016-3-18 22:17

楚香饭发表于 2015-9-30 07:35
由于浮点数有误差。所以 1.0/3.0 并不能精确的表示 1/3，例如可能是 0.333333334，极大可能是一个偶数分之 ...

麻烦问一下师兄，RRL=SQRT(UUAB**2+CH**2*(FF1**2-FF2**2))
RRL0=(RRL-UUAB)/(CH*(FF1+FF2))
其中UUAB为正，CH为正

作者: wutong-sky 时间: 2016-3-18 22:22

楚香饭发表于 2015-9-30 07:35
由于浮点数有误差。所以 1.0/3.0 并不能精确的表示 1/3，例如可能是 0.333333334，极大可能是一个偶数分之 ...

RRL=SQRT(UUAB**2+CH**2*(FF1**2-FF2**2))
RRL0=(RRL-UUAB)/(CH*(FF1+FF2))
师兄，麻烦问一下，我debug在这两句这出现了float invalid operation，怎么改？谢谢

作者: vvt 时间: 2016-3-18 22:44
请参考楼上的案例，依葫芦画瓢。

作者: wutong-sky 时间: 2016-3-19 16:17

vvt 发表于 2016-3-18 22:44
请参考楼上的案例，依葫芦画瓢。

RRL0=(RRL-UUAB)/(CH*(FF1+FF2))
麻烦问一下，这一句怎么改？

作者: vvt 时间: 2016-3-19 20:11
float invalid operation 是运行时错误

运行时错误，属于动态问题。往往不仅仅是代码本身错误导致的。而是数据不合适，或算法不严谨导致的。

比如

[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode

integer n

write(*,*) "有100个苹果分给小朋友，"

write(*,*) "请输入小朋友的个数"

read(*,*) n

write(*,*) "每个小朋友获得" , 100/n , "个苹果"

end

这个代码本身是没有错误的，但是当用户输入 0 时，就会出现 float invalid operation

这就是动态错误。建议你通过 debug 单步调试查找这类错误的原因。

如果你不熟悉 debug 单步调试，可以看这里的视频教程：http://v.fcode.cn/video-debugger.html

作者: 深流水静水流深 时间: 2016-4-3 07:20
学习了，呵呵

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)

Powered by Discuz! X3.2