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
!===================================================================================================================================

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