函数体构造对数值精度的过度反应
program mainimplicit none
real*8 :: Cd0,m,Vcruise,Rho,f1,f2,f3,f4,f5,f6,parent1,parent2,parent3,pi,funcval
! c=1.0d0
pi=4.0d0*datan(1.0d0)
Vcruise=150.0d0
m=1000.0d0
Cd0=0.02d0
Rho=1.212d0
funcval=0.0d0
parent1=11.34
parent2=0.67
parent3=32.83
f1=parent1/parent2!aspect ratio
f2=parent1*parent2!wing planform area
! f3=(((Wbp+4.936*f2*(f3**0.3)*9.81)*2)/(Rho*f2*(Vmin**2)))**2 !lift coefficient function (for minimum air speed)
f4=2.0d0*m*9.81d0/(Rho*f2*(Vcruise**2))!lift coefficient
if (parent3.le.30.0d0 .and. f1.gt.2.27d0 ) then
f5=1.78d0*(1.0d0-0.045d0*(f1**0.68d0))-0.64d0!Oswald Coefficient for sweep angle smaller than or equals to 30 Degree,f1 cannot be smaller than 2.27
elseif (parent3.gt.30.0d0 .and. f1.gt.10.0d0) then
f5=4.61d0*(1.0d0-0.045d0*(f1**0.68d0))*((dcos(parent3*pi/180.0d0))**0.15d0)-3.1d0!Oswald Coefficient for sweep angle larger than 30 Degree
else
f5=0.8d0
endif
f6=Cd0+(f4**2)/(pi*f1*f5)!drag coefficient
funcval=dble(Rho*(Vcruise**2)*f2*f6/2)
write(*,*) funcval
pause
end
以上函数实例，目前parent1，parent2，parent3给出的数值都为双精度小数点后两位有效数字，此时计算结果正确显示。
如果尝试人为或程式化直接产生parent1，parent2，parent3的双精度或单精度数值，则计算结果非正常显示。
此例子说明了在fortran编译器中，当函数或算式分母可能是非常小的实数时，参数的类型及计算误差会造成整体计算错误，所以应尽量回避这种情况。比如F=1/f。另外还可以人为添加约束条件。针对本例，主要问题集中于函数f5，对f5加入适当约束，比如f5小于0或一个给定值的时候，强制给定f5等于一个较大的数，可以有效解决计算失效的问题；或者可以将对应参数parent1，parent2，parent3的值强制截断为小数点后保留两位有效数字，命令如下:
parent(i)=(dint(parent(i)*100.0d0))/100.0d0
没看明白 函数计算的时候，参数的数据类型和精度会影响函数的计算结果，尤其是当作除法的时候，如果分子很大，分母很小，那么会出现计算结果不正确的情况。而此时编译器在编译执行过程中并不报错。我帖子里描述的很清楚了，你自己赋值黏贴代码，然后自己编译执行以下，之后自己在三个参数parent的数值后面随便加一些不是0的数字，按双精度长度加，之后重新编译执行，就会看到错误结果。加的越多，错误越离谱。 Felix804665 发表于 2014-3-5 01:41
函数计算的时候，参数的数据类型和精度会影响函数的计算结果，尤其是当作除法的时候，如果分子很大，分母很 ...
因为计算机字长的关系，相差悬殊的数字的算术运算是会有偏差的，但一般规律是采用的kind的精度越高计算结果应该越准确。 现在一般都定义双精度了。。。。
页:
[1]