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

[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