0

Please I need help with fixing the "unclassifiable statement error" I get when I try to call a function rtsafe, which finds the root of a nonlinear algebraic equation; using the Newton-Raphson and bisection method. The function takes in three inputs: x1 and x2, which are initial guesses in which the root lies between them, and xacc the third input is the tolerance used in the convergence test. When I compile and run the function alone, it doesn't give any error. However, when I call the function in the program, it gives the error:

   rtsafe(0.09,1.0,0.001)
      1
Error: Unclassifiable statement at (1)

kindly find the code below. Please note that I have attached to the code the FUNCTION which returns the given function and it's derivative at value x:

      program rootfinding

      implicit none
      real(4) :: rtsafe

      rtsafe(0.09,1.0,0.001)
      
      end program

      function rtsafe(x1,x2,xacc)
      implicit none
      integer,parameter :: maxit = 100
      real(4) :: x1,x2, xacc,
      real(4) :: rtsafe
      integer :: j
      real(4) :: df, dx, dxold, fh, fl, temp, xh, xl, f
      real(8):: 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

c     function for the function at value x
      function func(x)
      implicit none
      real(8):: func
      real(4) :: x
      func = x**2 -0.25
      end function

c     function for the derivative of the function at value x
      function funcd(x)
      implicit none
      real(4):: funcd
      real(4) :: x
      funcd= 2*x
      end function

I tried to check if there was any improper use of the 'if ' statement in the code but couldn't find any.

Tommy
  • 1
  • 2
  • `rtsafe` being a function, you have to use it inside an expression (e.g. `foo = rtsafe(0.09,1.0,0.001)`, or `write(*,*) rtsafe(0.09,1.0,0.001)`) – PierU Dec 26 '22 at 22:50

0 Answers0