FUNCTION zeroin(f, ax, bx, aerr, rerr) RESULT(fn_val) ! Code converted using TO_F90 by Alan Miller ! Date: 2003-07-14 Time: 12:32:54 !----------------------------------------------------------------------- ! FINDING A ZERO OF THE FUNCTION F(X) IN THE INTERVAL (AX,BX) ! ------------------------ ! INPUT... ! F FUNCTION SUBPROGRAM WHICH EVALUATES F(X) FOR ANY X IN THE ! CLOSED INTERVAL (AX,BX). IT IS ASSUMED THAT F IS CONTINUOUS, ! AND THAT F(AX) AND F(BX) HAVE DIFFERENT SIGNS. ! AX LEFT ENDPOINT OF THE INTERVAL ! BX RIGHT ENDPOINT OF THE INTERVAL ! AERR THE ABSOLUTE ERROR TOLERANCE TO BE SATISFIED ! RERR THE RELATIVE ERROR TOLERANCE TO BE SATISFIED ! OUTPUT... ! ABCISSA APPROXIMATING A ZERO OF F IN THE INTERVAL (AX,BX) !----------------------------------------------------------------------- ! ZEROIN IS A SLIGHTLY MODIFIED TRANSLATION OF THE ALGOL PROCEDURE ! ZERO GIVEN BY RICHARD BRENT IN ALGORITHMS FOR MINIMIZATION WITHOUT ! DERIVATIVES, PRENTICE-HALL, INC. (1973). !----------------------------------------------------------------------- IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) REAL (dp), INTENT(IN) :: ax REAL (dp), INTENT(IN) :: bx REAL (dp), INTENT(IN) :: aerr REAL (dp), INTENT(IN) :: rerr REAL (dp) :: fn_val ! EXTERNAL f INTERFACE FUNCTION f(x) RESULT(fn_val) IMPLICIT NONE INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(14, 60) REAL (dp), INTENT(IN) :: x REAL (dp) :: fn_val END FUNCTION f END INTERFACE REAL (dp) :: a, b, c, d, e, eps, fa, fb, fc, tol, xm, p, q, r, s, atol, rtol ! COMPUTE EPS, THE RELATIVE MACHINE PRECISION eps = EPSILON(0.0_dp) ! INITIALIZATION a = ax b = bx fa = f(a) fb = f(b) atol = 0.5 * aerr rtol = MAX(0.5_dp*rerr, 2.0_dp*eps) ! BEGIN STEP 10 c = a fc = fa d = b - a e = d 20 IF (ABS(fc) < ABS(fb)) THEN a = b b = c c = a fa = fb fb = fc fc = fa END IF ! CONVERGENCE TEST tol = rtol * MAX(ABS(b),ABS(c)) + atol xm = 0.5 * (c-b) IF (ABS(xm) > tol) THEN IF (fb /= 0.0) THEN ! IS BISECTION NECESSARY IF (ABS(e) >= tol) THEN IF (ABS(fa) > ABS(fb)) THEN ! IS QUADRATIC INTERPOLATION POSSIBLE IF (a == c) THEN ! LINEAR INTERPOLATION s = fb / fc p = (c-b) * s q = 1.0 - s ELSE ! INVERSE QUADRATIC INTERPOLATION q = fa / fc r = fb / fc s = fb / fa p = s * ((c-b)*q*(q-r)-(b-a)*(r-1.0)) q = (q-1.0) * (r-1.0) * (s-1.0) END IF ! ADJUST SIGNS IF (p > 0.0) q = -q p = ABS(p) ! IS INTERPOLATION ACCEPTABLE IF (2.0*p < (3.0*xm*q-ABS(tol*q))) THEN IF (p < ABS(0.5*e*q)) THEN e = d d = p / q GO TO 30 END IF END IF END IF END IF ! BISECTION d = xm e = d ! COMPLETE STEP 30 a = b fa = fb IF (ABS(d) > tol) b = b + d IF (ABS(d) <= tol) b = b + SIGN(tol,xm) fb = f(b) IF (fb*(fc/ABS(fc)) > 0.0) GO TO 10 GO TO 20 END IF END IF ! DONE fn_val = b RETURN END FUNCTION zeroin