NewtonRoot1D Function

public function NewtonRoot1D(tol, a, b, maxstep, xin, F0, fobj) result(xout)

Newton's iterative algorithm to find the root of function FR(x(:)) in the interval [a(:),b(:)], using d/dx(:)F(x)=0 and the derivative

Arguments

Type IntentOptional Attributes Name
real(kind=wp), intent(in) :: tol

abort tolerance
real(kind=wp), intent(in) :: a

search interval
real(kind=wp), intent(in) :: b

search interval
real(kind=wp), intent(in) :: maxstep

max|dx| allowed
real(kind=wp), intent(in) :: xin

initial guess
real(kind=wp), intent(in) :: F0

function to find root is FR(x)-F0
class(c_newton_Root1D) :: fobj

Return Value real(kind=wp)

on output =f(x)


Calls

proc~~newtonroot1d~~CallsGraph proc~newtonroot1d NewtonRoot1D FR FR proc~newtonroot1d->FR dFR dFR proc~newtonroot1d->dFR

Source Code

FUNCTION NewtonRoot1D(tol,a,b,maxstep,xin,F0,fobj) RESULT (xout)
  ! MODULES
  IMPLICIT NONE
  ! INPUT VARIABLES -------------------------!
  REAL(wp),INTENT(IN) :: tol    !! abort tolerance
  REAL(wp),INTENT(IN) :: a,b    !! search interval
  REAL(wp),INTENT(IN) :: maxstep !! max|dx| allowed
  REAL(wp),INTENT(IN) :: xin    !! initial guess
  REAL(wp),INTENT(IN) :: F0     !! function to find root is FR(x)-F0
  ! OUTPUT VARIABLES -------------------------!
  CLASS(c_newton_Root1D) :: fobj
  REAL(wp)            :: xout    !! on output =f(x)
  ! LOCAL VARIABLES --------------------------!
  INTEGER             :: iter,maxiter
  REAL(wp)            :: x,dx
  LOGICAL             :: converged
  LOGICAL             :: converged2
  ! CODE --------------------------------------------------------------------------------------------------------------------------!

  converged=.FALSE.
  x=xin
  maxiter=20
  DO iter=1,maxiter
    dx=-(fobj%FR(x)-F0)/fobj%dFR(x)
    dx = MAX(-(x-a),MIN(b-x,dx)) !respect bounds
    x = x+dx
    IF(ABS(dx).GT.maxstep) dx=dx/ABS(dx)*maxstep
    converged=(ABS(dx).LT.tol).AND.(x.GT.a).AND.(x.LT.b)
    IF(converged) EXIT
  END DO !iter
  IF(.NOT.converged) THEN
    !repeat with maxstep /10 and a little change in the initial condition
    x=MIN(b,MAX(a,xin+0.01_wp*(b-a)))
    maxiter=200
    DO iter=1,maxiter
      dx=-(fobj%FR(x)-F0)/fobj%dFR(x)
      dx = MAX(-(x-a),MIN(b-x,dx)) !respect bounds
      IF(ABS(dx).GT.maxstep) dx=dx/ABS(dx)*0.1_wp*maxstep
      x = x+dx
      converged2=(ABS(dx).LT.tol).AND.(x.GT.a).AND.(x.LT.b)
      IF(converged2) EXIT
    END DO !iter
    IF(converged2) THEN
      xout=x
      RETURN
    END IF
    WRITE(UNIT_stdout,*)'Newton abs(dx)<tol',ABS(dx),tol,ABS(dx).LT.tol
    WRITE(UNIT_stdout,*)'Newton x>a',x,a,(x.GT.a)
    WRITE(UNIT_stdout,*)'Newton x<b',x,b,(x.LT.b)
    WRITE(UNIT_stdout,*)'after iter',iter-1
    CALL abort(__STAMP__, &
              'NewtonRoot1D not converged')
  END IF
  xout=x

END FUNCTION NewtonRoot1D