关于二分法出现的奇怪问题求助

coldfield

新人求助，自己根据书上例子写了一段二分法计算的程序，有些数据可以计算出结果，有些却陷入死循环，这些数据相差不大并没有什么异常，在例子中1-23组的数据可以正常计算出，其它数据却不行。
求大神看看这段程序哪里有问题，万分感激！

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

!---------------------------------------------------------------------
!定义网格文件数据类型
!---------------------------------------------------------------------
Module Grid_data
  Implicit None
    Type Node
      Real(Kind=4) :: Radius, Length, Width, Waterbuild_up, WallShear
      Real(Kind=4) :: Valu
    End Type

  Real, Parameter :: PAI=3.1415926535 ,&
                                WV=1.155E-3          ,&
                                WD=999.0               ,&
                                N=17188.7              ,&
                                Omega=2.0*PAI*N/60.0  
                                 
  Real, Parameter :: RBLADE=36.0
End Module
!---------------------------------------------------------------------------
!定义方程
!---------------------------------------------------------------------------
Module WaterMass
  Implicit None
    Real(Kind=4) :: Y
    Real(Kind=4) :: C1, C2
    Real(Kind=4) :: V

  Contains
  Real Function F(X)
    Implicit None
    Real(Kind=4) :: X
    F=X+C1*C2*X**(2.0/3.0)-Y
    Return
  End Function
End Module
!------------------------------------------------------------------------------
!利用二分法求解方程
!-------------------------------------------------------------------------------
Module Numerical
Use WaterMass
  Implicit None
  Real, Parameter :: zero=1.0E-6

Contains
  Real Function Get_MK(P, Q, Func)
    Implicit None
    Real(Kind=4) :: P, Q  
    Real(Kind=4) :: R
    Real(Kind=4) :: FP
    Real(Kind=4) :: FQ
    Real(Kind=4) :: FR
    Real(Kind=4), External :: Func

    R=(P+Q)/2.0
    FR=Func(R)
    Do While(ABS(FR)>zero)   !-------------------------------------------------------------------------------------------------
      FP=Func(P)                       !经检测这里出现了奇怪的现象，在Do While语句之前并没有给FP、FQ赋值，FP、FQ、P、Q的值正常
      FQ=Func(Q)                      !经过Do While判断语句后，P和Q，FP和FQ的值莫名的相等了，使循环进入死循环
      If(FP*FR<0.0) Then           !但是，有些数据可以正常算出，有些相差不大的数据却陷入死循环，这点真的很郁闷⊙︿⊙
        Q=R                                   !--------------------------------------------------------------------------------------------------
        R=(P+Q)/2.0
      Else
        P=R
        R=(P+Q)/2.0
      End If
        FR=Func(R)
    End Do
    Get_MK=R
    Return
  End Function
End Module

Program Film_motion
  Use Grid_data
  Use Numerical
  Implicit None
  Integer, Parameter :: Fileid=100
  Integer, Parameter :: Nodes=295
  Character(len=80) :: Tempstr
  Type(Node) :: A(Nodes)
  Type(Node) :: B(Nodes)
  Type(Node) :: MK(Nodes)
  Type(Node) :: Valu
  Integer i, num, error

  Open(Fileid, File="build-up water_test.txt",status="old", iostat=error)
  If (error/=0) Then
    Write(*,*)"Open build-up water_test.txt fail."
    Stop
  End If

  Read(Fileid, "(A80)")Tempstr
  Valu=Node(0.0,0.0,0.0,0.0,0.0,0.0)
  !----------------------------------------------------------------------------------------------
  !读取文件数据，计算每个网格节点处A，B的值，A、B的值可以正常计算出
  !----------------------------------------------------------------------------------------------
  Do i=1, Nodes
    Read(Fileid,*)num, B(i)%Radius, B(i)%Length, B(i)%Width,&
                        B(i)%Waterbuild_up, B(i)%WallShear
    B(i)%Valu=B(i)%Length**(5.0/3.0)*B(i)%WallShear/((B(i)%Length*B(i)%Width*&
                     B(i)%Radius))**(2.0/3.0)
    A(i)%Valu=(WD/(2.0*WV))*(3.0*WV/(WD**2.0*Omega**2.0))**(2.0/3.0)
  End Do

  Write(*,"(2XA5,A15,4XA15,8XA15)")"Node","Radius", "B","A"
  Do i=1, Nodes
    Write(*,"(I5,3E22.10)") i, B(i)%Radius, B(i)%Valu , A(i)%Valu
  End Do

Pause
  !--------------------------------------------------------------------------------------------------
  !调用二分法函数，计算每个网格节点处MK的值，MK的计算出现异常，原因初步定为二分法函数的问题
  !--------------------------------------------------------------------------------------------------
  Do  i=1, Nodes
    C1=A(i)%Valu
    C2=B(i)%Valu
    Y=B(i)%Waterbuild_up
    V=0.0
    If(F(V)*F(Y)<0.0) Then
      MK(i)%Valu=Get_MK(V, Y, F)
    Else
      MK(i)%Valu=0.0
    End If
  End Do

  Write(*,"(1XA5,2XA10,4XA10,10XA10,10XA10)")"Node","Radius", "B","A", "MK"
  Do i=1, Nodes
    Write(*,"(I5,4E18.10)") i, B(i)%Radius, B(i)%Valu , A(i)%Valu, MK(i)%Valu
  End Do

  Stop
End Program

fcode

所有的数值解法，都有一定的适用范围。
二分法是有这个问题的。当 deltaX 小到低于浮点数有效位数（单精度为6-7位），还不足以引起大于 zero 的 deltaY 变化时，就会导致 P 和 Q 的值相同。
所以一般二分法在实际应用时，会增加一个终止条件，也就是判断 P 和 Q 的差值，如果小到一定的程度就终止，或循环到一定次数终止。

如果只用一个 FR 做为终止条件，是很危险的。

在很多迭代法中，都存在相同的问题。

coldfield

fcode 发表于 2015-10-6 21:00
所有的数值解法，都有一定的适用范围。
二分法是有这个问题的。当 deltaX 小到低于浮点数有效位数（单精度 ...

太感谢了，问题解决了，在ABS(FR)>zero后面加一个.AND. ABS（P-Q）>zero都可以算出来了，刚接触编程完全还没考虑到浮点数的精度的问题，非常感谢！
再问一下，如果把浮点数精度Kind值都改为8为双精度，提高浮点数有效位数或降低zero的值这样也可以算出结果吧？

fcode

提高浮点数精度不一定能解决你的问题（我没试过，但是这主要跟你的函数【的导数】有关）。

降低 zero 的要求应该可以。

coldfield

fcode 发表于 2015-10-6 21:42
提高浮点数精度不一定能解决你的问题（我没试过，但是这主要跟你的函数【的导数】有关）。

降低 zero 的要 ...

嗯，非常感谢，刚才试了将zero改为1.0E-1的量级也能算，算出的结果跟上个方法完全一样，这样的结果是不是不够精确？