函数体构造对数值精度的过度反应

Felix804665 · 发表于 2014-3-4 19:08:06

[Fortran] 纯文本查看 复制代码

program main

implicit 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

aliouying · 发表于 2014-3-4 23:47:50

没看明白

Felix804665 · 发表于 2014-3-5 01:41:46

函数计算的时候，参数的数据类型和精度会影响函数的计算结果，尤其是当作除法的时候，如果分子很大，分母很小，那么会出现计算结果不正确的情况。而此时编译器在编译执行过程中并不报错。我帖子里描述的很清楚了，你自己赋值黏贴代码，然后自己编译执行以下，之后自己在三个参数parent的数值后面随便加一些不是0的数字，按双精度长度加，之后重新编译执行，就会看到错误结果。加的越多，错误越离谱。