NewtonRoot2D Function

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

Newton's iterative algorithm to find the root of function [f1(x1,x2),f2(x1,x2)]=[0,0] in the interval a(i)<=x(i)<=b(i), using the Jacobian dfi/dxj, i=1,2, j=1,2, such that fi(x1,x2)=fi(x1_0,x2_0)+ [dfi/dx1,dfi/dx2].[dx1,dx2] in each step, we find dx1,dx2 st -[[dfi/dxj]] dxj =fi(x1_0,x2_0)

Arguments

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

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

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

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

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

initial guess
class(c_newton_Root2D) :: fobj

function to find root f1(x1,x2)=0,f2(x1,x2)=0 and derivatives d fi(x1,x2)/dxj

Return Value real(kind=wp), (2)

x1,x2 that have f1(x1,x2)=0 and f2(x1,x2)=0


Calls

proc~~newtonroot2d~~CallsGraph proc~newtonroot2d NewtonRoot2D FR FR proc~newtonroot2d->FR dFR dFR proc~newtonroot2d->dFR

Source Code

FUNCTION NewtonRoot2D(tol,a,b,maxstep,xin,fobj) RESULT (xout)
! MODULES
IMPLICIT NONE
!-----------------------------------------------------------------------------------------------------------------------------------
! INPUT VARIABLES
REAL(wp),INTENT(IN)    :: tol        !! abort tolerance
REAL(wp),INTENT(IN)    :: a(2),b(2)  !! search interval (2D)
REAL(wp),INTENT(IN) :: maxstep(2) !! max|dx| allowed
REAL(wp),INTENT(IN)    :: xin(2) !! initial guess
!-----------------------------------------------------------------------------------------------------------------------------------
! OUTPUT VARIABLES
CLASS(c_newton_Root2D) :: fobj !! function to find root f1(x1,x2)=0,f2(x1,x2)=0 and derivatives d fi(x1,x2)/dxj
REAL(wp)               :: xout(2) !! x1,x2 that have f1(x1,x2)=0 and f2(x1,x2)=0
!-----------------------------------------------------------------------------------------------------------------------------------
! LOCAL VARIABLES
INTEGER             :: iter,maxiter
REAL(wp)            :: x(2),dx(2)
REAL(wp)            :: det_Jac
REAL(wp)            :: F(2),Jac(2,2),JacInv(2,2)
LOGICAL             :: converged
!===================================================================================================================================
converged=.FALSE.
x=xin
maxiter=50
DO iter=1,maxiter
  Jac=fobj%dFR(x)
  det_Jac = Jac(1,1)*Jac(2,2)-Jac(1,2)*Jac(2,1)
  IF(det_Jac.LT.1.0E-12) CALL abort(__STAMP__,&
                                    'det Jacobian<=0 in NewtonRoot2d')
  JacInv(1,1)= Jac(2,2)
  JacInv(1,2)=-Jac(1,2)
  JacInv(2,1)=-Jac(2,1)
  JacInv(2,2)= Jac(1,1)
  JacInv=JacInv/det_Jac
  F=fobj%FR(x)
  dx=-MATMUL(JacInv,F)
  dx = MAX(-(x-a),MIN(b-x,dx)) !respect bounds
  IF(ABS(dx(1)).GT.maxstep(1)) dx(1)=dx(1)/ABS(dx(1))*maxstep(1)
  IF(ABS(dx(2)).GT.maxstep(2)) dx(2)=dx(2)/ABS(dx(2))*maxstep(2)

  x = x+dx
  converged=(SQRT(SUM(dx*dx)).LT.tol).AND.ALL(x.GT.a).AND.ALL(x.LT.b)
  IF(converged) EXIT
END DO !iter
xout=x
IF(.NOT.converged) THEN
  WRITE(UNIT_stdout,*)'Newton abs(dx)<tol',ABS(dx),tol,'F(x)',fobj%FR(xout)
  CALL abort(__STAMP__,&
             'NewtonRoot2D not converged')
END IF

END FUNCTION NewtonRoot2D