调用二分法子程序解方程出错，求助 - Fortran Coder

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

  Implicit none

  Integer :: i

  Real :: Maact1,Maact2

  Real,dimension(3) ::fact

  fact(1)=26.

  fact(2)=30.

  fact(3)=32.

  Do i=1,3                                      !输出的结果Maact1和Maact2结果一致，都为1.985840

    Call equation(fact(i),Maact1)               !但是把9-12行注释掉的话，

    Write(*,*) Maact1                           !输出的Maact2结果为2.1938，是正确的

  End Do

  Call equation(31.57,Maact2)

  Write(*,*) Maact2

  End Program test



  Subroutine equation(f,Ma)              !已知函数值，根据普朗特-迈耶函数求马赫数,1<=Ma<=3

  Implicit none

  Real,intent(in) :: f                   !输入普朗特-迈耶函数值

  Real,intent(out) :: Ma                 !输出马赫数

  Real :: eps=0.001                      !计算精度

  Real :: Ma1=1.0                        !区间下限

  Real :: Ma2=3.0                        !区间上限

  Real :: f_Ma1,f_Ma2,f_Ma               !对应函数值

  Real :: gamma=1.4

  Real,External :: P_Ma                  !普朗特迈耶函数

  !二分法解方程

  Ma=(Ma1+Ma2)/2

  f_Ma=P_Ma(Ma)-f

  f_Ma1=P_Ma(Ma1)-f

  f_Ma2=P_Ma(Ma2)-f

  Do while((abs(f_ma)>eps).and.((Ma2-Ma1)>eps))

    If(f_Ma*f_Ma1<0.) then

      Ma2=Ma

      Ma=(Ma1+Ma2)/2

      f_Ma2=f_Ma

    Else

      Ma1=Ma

      Ma=(Ma1+Ma2)/2

      f_Ma1=f_Ma

    End If

    f_Ma=P_Ma(Ma)-f

  End Do

  

  !f_Ma1=P_Ma(Ma1)-f

  !f_Ma2=P_Ma(Ma2)-f

  !Do

  !  Ma=(Ma1+Ma2)/2

  !  f_Ma=P_Ma(Ma)-f

  !  If ((abs(f_Ma)<=eps).or.abs(Ma1-Ma2)<=eps) exit

  !  If (f_Ma1*f_Ma<0) then

  !    Ma2=Ma

  !    f_Ma2=f_Ma

  !  Else

  !    Ma1=Ma

  !    f_Ma1=f_Ma

  !  End if

  !End do

  End subroutine equation



  Real Function P_Ma(Ma)

  Implicit None

  Real,Intent(in) :: Ma

  Real :: gamma=1.4

  P_Ma=sqrt((gamma+1)/(gamma-1))*atand(sqrt((gamma-1)/(gamma+1)*(Ma**2-1)))-atand(sqrt(Ma**2-1))

  End FUnction P_Ma