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