0

Please, I have this code which computes the root of a nonlinear algebraic equation. When I compile the code for the written function(rtsafe) alone, it works, but when I call the function rtsafe using y = rtsafe(0.01000000000,10.0000000000, 0.00000001) it complains about the argument input for 'x1', giving the error message:


      y= rtsafe(0.01000000000,10.0000000000,0.00000001)
                1
Warning: Type mismatch in argument 'x1' at (1); passed REAL(4) to REAL(8)

I don't understand how my 'X1 =0.01000000000' is REAL(4) because it has more than 8 decimal digits. kindly find below the code:


      program rootfinding
      real(8) :: rtsafe,y

      real(8) ::testroot

      y= rtsafe(0.01000000000,10.0000000000,0.00000001)
      
      testroot =(y**2 +0.25)**0.5-1
      print *, testroot,y

      end program

      function rtsafe(x1,x2,xacc)
      implicit none
      integer,parameter :: maxit = 200
      real(8) :: x1,x2
      real(8) :: rtsafe
      real(4) ::  xacc
      integer :: j
      real(4) :: df, dx, dxold, fh, fl, temp, xh, xl, f
      real(4):: func
      real(4):: funcd

      fl =  func(x1)
      fh = func(x2)
      if ((fl > 0 .and. fh > 0) .or. (fl < 0 .and. fh < 0)) then
         print *, 'root must be bracketed'
         return
      end if

      if (fl == 0) then
        rtsafe = x1
        return
      else if (fh == 0) then
        rtsafe = x2
        return

      c    orient the search so that f(x1) < 0

      else if (fl < 0) then
         xl = x1
         xh = x2
      else
         xh = x1
         xl = x2
      end if

      c    initialize  the guess for root,
         rtsafe = 0.5 * (x1+x2)

      c    the step size before last
         dxold = abs(x2 - x1)

      c    and the last step
         dx = dxold

      df = funcd(rtsafe)
      f = func(rtsafe)

      c    loop over allowed iterations

      do j = 1, maxit

      c    bisect if newton is out of range or not decreasing fast enough
          if (((rtsafe - xh)*df - f)*((rtsafe - xl)*df - f) > 0 .or.    &
     &           abs(2*f) > abs (dxold*df)) then
              dxold = dx
              dx = 0.5*(xh - xl)
              rtsafe = xl  + dx
              if (xl == rtsafe) return
          else
              dxold = dx
              dx = f/df
              temp = rtsafe
              rtsafe = rtsafe - dx
              if (temp == rtsafe ) return
          end if

      c     convergence criterion
         if (abs(dx) < xacc) return

         f = func(rtsafe)
         if (f < 0) then
              xl = rtsafe
         else
              xh = rtsafe
         end if

      end do
      print *, 'rtsafe exceeding maximum iterations'
      return
      end function rtsafe

      function func(x)
      implicit none
      real(4):: func
      real(8) :: x
      func = (x**2+ 0.25)**0.5-1.0
      end function

      function funcd(x)
      implicit none
      real(4):: funcd
      real(8) ::x
      funcd= x*(x**2 +0.25)**(-0.5)
      end function
Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
Tommy
  • 1
  • 2
  • 1
    A number like `0.001234567890123...` is **always** the default real kind, no matter how many digits it has. Please read more about how kinds in Fortran work and do not use the numbers 4 and 8, they are not portable https://stackoverflow.com/questions/838310/fortran-90-kind-parameter Also, always use tag [tag:fortran] for Fortran questions. – Vladimir F Героям слава Jan 13 '23 at 05:52
  • A similar problem https://stackoverflow.com/questions/24405760/fortran-passing-numbers-to-subroutines Here the same thing with integers https://stackoverflow.com/questions/31162311/type-of-a-hardcoded-argument – Vladimir F Героям слава Jan 13 '23 at 05:55
  • This one is really informative https://stackoverflow.com/questions/49183219/is-it-necessary-to-append-kind-to-number-literals-in-modern-fortran And if you work with numbers with many sigits, **definitely** consider problems like https://stackoverflow.com/questions/16370206/numerical-precision-in-fortran-95 https://stackoverflow.com/questions/6146005/is-there-a-better-double-precision-assignment-in-fortran-90?rq=1 – Vladimir F Героям слава Jan 13 '23 at 05:57

0 Answers0