Fortran Coder

标题: Fortran计算结果浮点数异常 [打印本页]

作者: 向着光亮 时间: 2018-1-22 20:13
标题: Fortran计算结果浮点数异常
如图。按照文献里的公式计算天文日照，只是纬度值稍微变了一点就计算有误，我把正确的和错误的计算的放上面了，不知道怎么回事，改了定义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, 下载次数: 455)

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

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

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

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

作者: pasuka 时间: 2018-1-22 22:16
本帖最后由 pasuka 于 2018-1-22 22:18 编辑

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

作者: vvt 时间: 2018-1-23 09:42
代码复制粘贴更好，不要截图

作者: 向着光亮 时间: 2018-1-23 09:44

vvt 发表于 2018-1-23 09:42
代码复制粘贴更好，不要截图

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

作者: 向着光亮 时间: 2018-1-23 09:56

pasuka 发表于 2018-1-22 22:16
请将公式的系数，譬如1.00011表示成单精度1.00011E0，双精度1.00011D0，180改成单精度1.8E2，双精度1.8D2
...

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

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

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

作者: vvt 时间: 2018-1-23 10:08

向着光亮发表于 2018-1-23 09:44
为什么呢，这样也能看清楚啊

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

作者: 向着光亮 时间: 2018-1-23 10:12

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

作者: vvt 时间: 2018-1-23 10:40
lat=42.4025*PI/180.0
纬度应该是弧度，而非角度。

作者: pasuka 时间: 2018-1-23 11:21

向着光亮发表于 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)?
师傅领进门，修行在自身

作者: 向着光亮 时间: 2018-1-23 11:22

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, 下载次数: 390)

作者: 向着光亮 时间: 2018-1-23 11:25

vvt 发表于 2018-1-23 10:40
lat=42.4025*PI/180.0
纬度应该是弧度，而非角度。

谢谢你，我明白了原理，非常感谢，已经算出来了

作者: 向着光亮 时间: 2018-1-23 11:27

pasuka 发表于 2018-1-23 11:21
依旧是问题代码：
1、明明打开同一个文件，open语句为何放在do循环内？
2、有open无close；

你说的对，这样也能算出来，把open放出来和在里面都是可以的，但是习惯很重要，所以我决定拿出来了，，至于这个单双精度浮点问题刚开始接受，需要慢慢适应

作者: 向着光亮 时间: 2018-1-23 11:36

pasuka 发表于 2018-1-23 11:21
依旧是问题代码：
1、明明打开同一个文件，open语句为何放在do循环内？
2、有open无close；

你好，听了你的建议我觉得确实应该严谨的写才能避免出错，我还想问问你这个计算的时候所有的数字都得写的浮点型的吗，定义变量也要这样吗

作者: 向着光亮 时间: 2018-1-23 12:40

vvt 发表于 2018-1-23 10:40
lat=42.4025*PI/180.0
纬度应该是弧度，而非角度。

我刚才又算了一遍，还是不行，大体思路我明白了，主要是角度和弧度的互换，太容易出错了，你刚才说的纬度需要换成弧度我觉得不对，我重新编了代码，后面都表明了单位，可是还是浮点数异常或者NAN，实在不知为何，计算了一上午，求解

[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*pi)*(n-1)/366!年角:弧度



 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)=ws(n)*(2/15)*(pi/180)!日长：弧度

end do

!write(*,fmt="(365f8.2)")delta(n)

open(unit=7,file="D:\Fortran读写文件\计算云量\tiwen\2012tianwen")

write(*,*)sl*10



stop

end

QQ截图20180123124009.png (50.06 KB, 下载次数: 358)

作者: fcode 时间: 2018-1-23 12:58
tan 就是弧度制的，不是你想什么制，就什么制。

作者: pasuka 时间: 2018-1-23 13:21

向着光亮发表于 2018-1-23 12:40
我刚才又算了一遍，还是不行，大体思路我明白了，主要是角度和弧度的互换，太容易出错了，你刚才说的纬度 ...

坚持角度制的话，请查找手册，譬如：
https://gcc.gnu.org/onlinedocs/g ... an/ASIND.html#ASIND

9.20 ASIND — Arcsine function, degrees
Description:
ASIND(X) computes the arcsine of its X in degrees (inverse of SIND(X)).

This function is for compatibility only and should be avoided in favor of standard constructs wherever possible.

Standard:
GNU Extension, enabled with -fdec-math.

Class:
Elemental function

Syntax:
RESULT = ASIND(X)

Arguments:
X The type shall be either REAL and a magnitude that is less than or equal to one - or be COMPLEX.
Return value:
The return value is of the same type and kind as X. The real part of the result is in degrees and lies in the range -90 \leq \Re \asin(x) \leq 90.

Example:
program test_asind
real(8) :: x = 0.866_8
x = asind(x)
end program test_asind
Specific names:
Name Argument Return type Standard
ASIND(X) REAL(4) X REAL(4) GNU Extension
DASIND(X) REAL(8) X REAL(8) GNU Extension
See also:
Inverse function: SIND Radians function: ASIN

作者: 向着光亮 时间: 2018-1-23 21:18

fcode 发表于 2018-1-23 12:58
tan 就是弧度制的，不是你想什么制，就什么制。

tan计算的不就是一个数值吗？为何是弧度制？

作者: pasuka 时间: 2018-1-23 21:49

向着光亮发表于 2018-1-23 21:18
tan计算的不就是一个数值吗？为何是弧度制？

请先回顾计算方法课程的相关内容：三角函数值
传送门：
http://blog.csdn.net/eric88/article/details/17031311

作者: vvt 时间: 2018-1-23 22:59

向着光亮发表于 2018-1-23 21:18
tan计算的不就是一个数值吗？为何是弧度制？

tan 的定义域是以弧度制为单位，而不是值域。
所以，tan(lat) 要求 lat 是以弧度制为单位。其结果（值域），是个比值，无单位。

作者: 向着光亮 时间: 2018-1-24 18:25

vvt 发表于 2018-1-23 22:59
tan 的定义域是以弧度制为单位，而不是值域。
所以，tan(lat) 要求 lat 是以弧度制为单位。其结果（值域 ...

好的，明白了

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