Question

【急】fortran 运行结果出现NaN

一道fortran程序题，求解一个函数的解，但是结果一直是nan，但是从数学上看这个程序编得没有问题，换一种写法后结果可以出来，但是我很想知道为什么这样写会错（编译没有错误），望高手解答！源程序如下：
function f(x) result (s)
implicit none
real::x,un1,un2,un,s,eps=1e-10
integer::n
un=1.0
un2=1.0
un1=1.0
s=0.5*cos(3*x)
n=2
do while(abs(un)>eps)
un1=un1*(2*n-1)
un2=un2*(2*n)
un=un1/un2*cos((2*n+1)*x)
s=s+un
n=n+1
enddo
s=s+cos(x)
end function

program main
implicit none
real::x,f
integer::i
character(len=60)::fmt="('x=',f5.3,'f(x)=',f8.4)"
open(1,file='myf2.dat')
open(2,file='myf2.out')

do i=1,7
read(1,*)x

write(*,fmt)x,f(x)
write(2,fmt)x,f(x)
enddo
end

跪谢
我不明白一楼说的会出现无穷。
我在循环处设了范围界定，un的大小被限制住了，这样n也不可能无穷大啊，n不是无穷大，那un1=un1*(2*n-1)也就不应该会无穷啊

Answer 1

原因在于楼主要算的级数常常要在30项开外才收敛，因此
un1=3*5*7*...*(2*30-1)*...
un2=4*6*8*...*(2*30)*...
超出单精度浮点数的范围3.40E38

解决方法是去掉un1，un2，直接计算un。
! un1=un1*(2*n-1)
! un2=un2*(2*n)
un=un*(2*n-1)/(2*n)*cos((2*n+1)*x)

我试了，即使un1 un2改成双精度也不行。比如x=1.0时，需要152项才收敛。而算到第150项时，un1恰好达到双精度浮点数的上界1E308。

Answer 2

我不清楚你myf2.dat文件里的数据时什么，但
un1=un1*(2*n-1)
un2=un2*(2*n)
这样算会出现无穷。

如果你那do while的循环计算30次，那么
un2=2^30*30!
这在计算机里是没法表示的。
你范围的界定是对un的绝对值进行界定，你可以试一下将eps值改大（如1e-2)就行了。