RECURSIVE SUBROUTINE SolveRegulaFalsi(Eps, MaxIte, Flag, XRes, f, X_0, X_1, Par)
! SUBROUTINE INFORMATION:
! AUTHOR Michael Wetter
! DATE WRITTEN March 1999
! MODIFIED Fred Buhl November 2000, R. Raustad October 2006 - made subroutine RECURSIVE
! RE-ENGINEERED na
! PURPOSE OF THIS SUBROUTINE:
! Find the value of x between x0 and x1 such that f(x,[,Par])
! is equal to zero.
! METHODOLOGY EMPLOYED:
! Uses the Regula Falsi (false position) method (similar to secant method)
! REFERENCES:
! See Press et al., Numerical Recipes in Fortran, Cambridge University Press,
! 2nd edition, 1992. Page 347 ff.
! USE STATEMENTS:
! na
IMPLICIT NONE ! Enforce explicit typing of all variables in this routine
! SUBROUTINE ARGUMENT DEFINITIONS:
REAL(r64), INTENT(IN) :: Eps ! required absolute accuracy
INTEGER, INTENT(IN) :: MaxIte ! maximum number of allowed iterations
INTEGER, INTENT(OUT) :: Flag ! integer storing exit status
! = -2: f(x0) and f(x1) have the same sign
! = -1: no convergence
! > 0: number of iterations performed
REAL(r64), INTENT(OUT) :: XRes ! value of x that solves f(x [,Par]) = 0
REAL(r64), INTENT(IN) :: X_0 ! 1st bound of interval that contains the solution
REAL(r64), INTENT(IN) :: X_1 ! 2nd bound of interval that contains the solution
REAL(r64), DIMENSION(:), INTENT(IN), OPTIONAL :: Par ! array with additional parameters used for function evaluation
! optional
! SUBROUTINE PARAMETER DEFINITIONS:
REAL(r64), PARAMETER :: SMALL = 1.d-10
! INTERFACE BLOCK SPECIFICATIONS
INTERFACE ! Interface to function to be solved for zero: f(X, Par) = 0
FUNCTION f(X, Par) RESULT (Y)
USE DataPrecisionGlobals
REAL(r64), INTENT(IN) :: X
REAL(r64), INTENT(IN), DIMENSION(:), OPTIONAL :: Par
REAL(r64) :: Y
END FUNCTION
END INTERFACE
! DERIVED TYPE DEFINITIONS
! na
! SUBROUTINE LOCAL VARIABLE DECLARATIONS:
REAL(r64) :: X0 ! present 1st bound
REAL(r64) :: X1 ! present 2nd bound
REAL(r64) :: XTemp ! new estimate
REAL(r64) :: Y0 ! f at X0
REAL(r64) :: Y1 ! f at X1
REAL(r64) :: YTemp ! f at XTemp
REAL(r64) :: DY ! DY = Y0 - Y1
LOGICAL :: Conv ! flag, true if convergence is achieved
LOGICAL :: StopMaxIte ! stop due to exceeding of maximum # of iterations
LOGICAL :: Cont ! flag, if true, continue searching
INTEGER :: NIte ! number of interations
X0 = X_0
X1 = X_1
Conv = .FALSE.
StopMaxIte = .FALSE.
Cont = .TRUE.
NIte = 0
IF (PRESENT(Par)) THEN
Y0 = f(X0, Par)
Y1 = f(X1, Par)
ELSE
Y0 = f(X0)
Y1 = f(X1)
END IF
! check initial values
IF ( Y0*Y1 > 0 ) THEN
Flag = -2
XRes = X0
RETURN
END IF
DO WHILE (Cont)
DY = Y0 - Y1
IF (ABS(DY) < SMALL) DY = SMALL
! new estimation
XTemp = (Y0 * X1 - Y1 * X0 ) / DY
IF (PRESENT(Par)) THEN
YTemp = f(XTemp, Par)
ELSE
YTemp = f(XTemp)
END IF
NIte = NIte + 1
! check convergence
IF (ABS(YTemp) < Eps) Conv = .TRUE.
IF (NIte > MaxIte) StopMaxIte = .TRUE.
IF ((.NOT.Conv).AND.(.NOT.StopMaxIte)) THEN
Cont = .TRUE.
ELSE
Cont = .FALSE.
END IF
IF (Cont) THEN
! reassign values (only if further iteration required)
IF ( Y0 < 0.0d0 ) THEN
IF ( YTemp < 0.0d0 ) THEN
X0 = XTemp
Y0 = YTemp
ELSE
X1 = XTemp
Y1 = YTemp
END IF
ELSE
IF ( YTemp < 0.0d0 ) THEN
X1 = XTemp
Y1 = YTemp
ELSE
X0 = XTemp
Y0 = YTemp
END IF
END IF ! ( Y0 < 0 )
END IF ! (Cont)
END DO ! Cont
IF (Conv) THEN
Flag = NIte
ELSE
Flag = -1
END IF
XRes = XTemp
RETURN
END SUBROUTINE SolveRegulaFalsi