Fortran Coder

查看: 31412|回复: 19
打印 上一主题 下一主题

[数值问题] Fortran计算结果浮点数异常

[复制链接]

50

帖子

11

主题

0

精华

熟手

F 币
239 元
贡献
151 点
跳转到指定楼层
楼主
发表于 2018-1-22 20:13:18 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式
如图。按照文献里的公式计算天文日照,只是纬度值稍微变了一点就计算有误,我把正确的和错误的计算的放上面了,不知道怎么回事,改了定义real*8:: delta(365),y(365),ws(365),sl(365)为单双精度都试了,结果还是Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG,然后把纬度改为41.25它又有结果了,实在不知为何,忘有人指导!先行谢过

QQ截图20180122195648.png (21.91 KB, 下载次数: 534)

QQ截图20180122195648.png

QQ截图20180122195707.png (24.02 KB, 下载次数: 570)

QQ截图20180122195707.png

QQ截图20180122195752.png (22.98 KB, 下载次数: 553)

QQ截图20180122195752.png

QQ截图20180122195802.png (42.76 KB, 下载次数: 539)

QQ截图20180122195802.png

QQ截图20180122200950.png (15 KB, 下载次数: 550)

QQ截图20180122200950.png
分享到:  微信微信
收藏收藏 点赞点赞 点踩点踩

490

帖子

4

主题

0

精华

大宗师

F 币
3298 元
贡献
1948 点

水王勋章元老勋章热心勋章

沙发
发表于 2018-1-22 22:16:42 | 只看该作者
本帖最后由 pasuka 于 2018-1-22 22:18 编辑

请将公式的系数,譬如1.00011表示成单精度1.00011E0,双精度1.00011D0,180改成单精度1.8E2,双精度1.8D2
年角代码中的(n-1)改写成real(n-1)

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
板凳
发表于 2018-1-23 09:42:30 | 只看该作者
代码复制粘贴更好,不要截图

50

帖子

11

主题

0

精华

熟手

F 币
239 元
贡献
151 点
地板
 楼主| 发表于 2018-1-23 09:44:36 | 只看该作者
vvt 发表于 2018-1-23 09:42
代码复制粘贴更好,不要截图

为什么呢,这样也能看清楚啊

50

帖子

11

主题

0

精华

熟手

F 币
239 元
贡献
151 点
5#
 楼主| 发表于 2018-1-23 09:56:54 | 只看该作者
pasuka 发表于 2018-1-22 22:16
请将公式的系数,譬如1.00011表示成单精度1.00011E0,双精度1.00011D0,180改成单精度1.8E2,双精度1.8D2
...

你好,按照你的方法改了之后可以计算,但是结果是空的,这是怎么回事

1.png (20.63 KB, 下载次数: 529)

1.png

2.png (2.86 KB, 下载次数: 464)

2.png

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
6#
发表于 2018-1-23 10:08:58 | 只看该作者
向着光亮 发表于 2018-1-23 09:44
为什么呢,这样也能看清楚啊

因为截图不能复制粘贴。别人懒得给你测试。
反正我没耐心照着敲一遍你的代码。

50

帖子

11

主题

0

精华

熟手

F 币
239 元
贡献
151 点
7#
 楼主| 发表于 2018-1-23 10:12:30 | 只看该作者
vvt 发表于 2018-1-23 10:08
因为截图不能复制粘贴。别人懒得给你测试。
反正我没耐心照着敲一遍你的代码。 ...

号的,我复制给你
[Fortran] 纯文本查看 复制代码
program ex0109
real*8:: delta(366),y(366),ws(366),sl(366)
integer n
real*8,parameter::pi=3.1415926
real*8 lat
lat=42.4025

do n=1,366
   y(n)=2.0*pi*real(n-1)/366.0!年角

 !delta(n)=0.006918-0.399912*cos(y(n))+0.070257*sin(y(n))-0.006758*cos(2*y(n))&
!+0.000907*sin(2*y(n))-0.002697*cos(3*y(n))+0.00148*sin(3*y(n))!太阳赤纬
delta(n)=6.918E-3-3.99912E-1*cos(y(n))+7.0257E-2*sin(y(n))-6.758E-3*cos(2*y(n))&
+9.07E-4*sin(2*y(n))-2.697E-3*cos(3*y(n))+1.48E-3*sin(3*y(n))

ws(n)=acos(-tan(lat)*tan(delta(n)))!时角(角度)
sl(n)=(1.8E2/pi)*ws(n)*(2.0E0/1.5E1)!日长
!write(*,fmt="(365f8.2)")delta(n)
open(unit=7,file="D:\Fortran读写文件\计算云量\tiwen\2012tianwen")
write(*,*)sl(n)
end do 

stop
end

954

帖子

0

主题

0

精华

大师

F 币
184 元
贡献
75 点

规矩勋章元老勋章新人勋章水王勋章热心勋章

QQ
8#
发表于 2018-1-23 10:40:15 | 只看该作者
lat=42.4025*PI/180.0  
纬度应该是弧度,而非角度。

490

帖子

4

主题

0

精华

大宗师

F 币
3298 元
贡献
1948 点

水王勋章元老勋章热心勋章

9#
发表于 2018-1-23 11:21:38 | 只看该作者
向着光亮 发表于 2018-1-23 10:12
号的,我复制给你
[mw_shl_code=fortran,true]program ex0109
real*8:: delta(366),y(366),ws(366),sl(366 ...

依旧是问题代码:
1、明明打开同一个文件,open语句为何放在do循环内?
2、有open无close;
3、为啥不能举一反三把3*y(n)这类表述也改成3.0E0*y(n)?
师傅领进门,修行在自身

50

帖子

11

主题

0

精华

熟手

F 币
239 元
贡献
151 点
10#
 楼主| 发表于 2018-1-23 11:22:04 | 只看该作者
vvt 发表于 2018-1-23 10:40
lat=42.4025*PI/180.0  
纬度应该是弧度,而非角度。

好像有道理,但是我图片中放的应该是计算的角度啊,你看看ws(n)=acos(-tan(lat)*tan(delta(n)))!时角

33.png (2.4 KB, 下载次数: 455)

33.png
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

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

GMT+8, 2024-12-27 09:37

Powered by Tencent X3.4

© 2013-2024 Tencent

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