Fortran Coder

查看: 17569|回复: 10
打印 上一主题 下一主题

[通用算法] 负数求三次方根的NaN问题,多谢啦!

[复制链接]

2

帖子

1

主题

0

精华

新人

精神病

F 币
13 元
贡献
6 点
跳转到指定楼层
楼主
发表于 2015-9-30 01:42:06 | 显示全部楼层 |只看大图 回帖奖励 |倒序浏览 |阅读模式
不说多的 直接上源代码
[Fortran] 纯文本查看 复制代码
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
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

2

帖子

1

主题

0

精华

新人

精神病

F 币
13 元
贡献
6 点
沙发
 楼主| 发表于 2015-9-30 08:30:05 | 显示全部楼层
楚香饭 发表于 2015-9-30 07:35
由于浮点数有误差。所以 1.0/3.0 并不能精确的表示 1/3,例如可能是 0.333333334,极大可能是一个偶数分之 ...

果然编译成功了 谢谢大神!
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-5-16 06:32

Powered by Tencent X3.4

© 2013-2024 Tencent

快速回复 返回顶部 返回列表