NewtonRoot1D_FdF Function

public function NewtonRoot1D_FdF(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 on input
real(kind=wp), intent(in) :: F0

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

function to find root f(x) & derivative d/dx f(x) as FRdFR method

Return Value real(kind=wp)

output x for f(x)=0


Calls

proc~~newtonroot1d_fdf~~CallsGraph proc~newtonroot1d_fdf NewtonRoot1D_FdF FRdFR FRdFR proc~newtonroot1d_fdf->FRdFR

Source Code

FUNCTION NewtonRoot1D_FdF(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 on input
REAL(wp),INTENT(IN) :: F0      !! function to find root is FR(x)-F0
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
CLASS(c_newton_Root1D_FdF) :: fobj !! function to find root f(x) & derivative d/dx f(x) as FRdFR method
REAL(wp)            :: xout    !! output x for f(x)=0
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER             :: iter,maxiter
REAL(wp)            :: x,dx
REAL(wp)            :: FRdFRx(2) !1: FR(x), 2: dFR(x)
LOGICAL             :: converged
LOGICAL             :: converged2
!===================================================================================================================================
converged=.FALSE.
x=xin
maxiter=20
DO iter=1,maxiter
  FRdFRx=fobj%FRdFR(x)
  dx=-(FRdFRx(1)-F0)/FRdFRx(2)
  dx = MAX(-(x-a),MIN(b-x,dx)) !respect bounds
  IF(ABS(dx).GT.maxstep) dx=dx/ABS(dx)*maxstep
  x = x+dx
  converged=(ABS(dx).LT.tol).AND.(x.GE.a).AND.(x.LE.b)
  IF(converged) EXIT
END DO !iter
IF(.NOT.converged) THEN
  !repeat with maxstep /10 and a little change in the initial condition
  converged2=.FALSE.
  x=MIN(b,MAX(a,xin+0.01_wp*(b-a)))
  maxiter=200
  DO iter=1,maxiter
    FRdFRx=fobj%FRdFR(x)
    dx=-(FRdFRx(1)-F0)/FRdFRx(2)
    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.GE.a).AND.(x.LE.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_FdF not converged')
END IF
xout=x

END FUNCTION NewtonRoot1D_FdF