SUBROUTINE FZERO(F,B,C,R,RE,AE,IFLAG) C***BEGIN PROLOGUE FZERO C***DATE WRITTEN 700901 (YYMMDD) C***REVISION DATE 860411 (YYMMDD) C***CATEGORY NO. F1B C***KEYWORDS BISECTION,NONLINEAR,ROOTS,ZEROS C***AUTHOR SHAMPINE,L.F.,SNLA C WATTS,H.A.,SNLA C***PURPOSE FZERO searches for a zero of a function F(X) in a given C interval (B,C). It is designed primarily for problems C where F(B) and F(C) have opposite signs. C***DESCRIPTION C C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C Based on a method by T J Dekker C written by L F Shampine and H A Watts C C FZERO searches for a zero of a function F(X) between C the given values B and C until the width of the interval C (B,C) has collapsed to within a tolerance specified by C the stopping criterion, ABS(B-C) .LE. 2.*(RW*ABS(B)+AE). C The method used is an efficient combination of bisection C and the secant rule. C C Description Of Arguments C C F,B,C,R,RE and AE are input parameters C B,C and IFLAG are output parameters (flagged by an * below) C C F - Name of the real valued external function. This name C must be in an EXTERNAL statement in the calling C program. F must be a function of one real argument. C C *B - One end of the interval (B,C). The value returned for C B usually is the better approximation to a zero of F. C C *C - The other end of the interval (B,C) C C R - A (better) guess of a zero of F which could help in C speeding up convergence. If F(B) and F(R) have C opposite signs, a root will be found in the interval C (B,R); if not, but F(R) and F(C) have opposite C signs, a root will be found in the interval (R,C); C otherwise, the interval (B,C) will be searched for a C possible root. When no better guess is known, it is C recommended that r be set to B or C; because if R is C not interior to the interval (B,C), it will be ignored. C C RE - Relative error used for RW in the stopping criterion. C If the requested RE is less than machine precision, C then RW is set to approximately machine precision. C C AE - Absolute error used in the stopping criterion. If the C given interval (B,C) contains the origin, then a C nonzero value should be chosen for AE. C C *IFLAG - A status code. User must check IFLAG after each call. C Control returns to the user from FZERO in all cases. C C 1 B is within the requested tolerance of a zero. C The interval (B,C) collapsed to the requested C tolerance, the function changes sign in (B,C), and C F(X) decreased in magnitude as (B,C) collapsed. C C 2 F(B) = 0. However, the interval (B,C) may not have C collapsed to the requested tolerance. C C 3 B may be near a singular point of F(X). C The interval (B,C) collapsed to the requested tol- C erance and the function changes sign in (B,C), but C F(X) increased in magnitude as (B,C) collapsed,i.e. C abs(F(B out)) .GT. max(abs(F(B in)),abs(F(C in))) C C 4 No change in sign of F(X) was found although the C interval (B,C) collapsed to the requested tolerance. C The user must examine this case and decide whether C B is near a local minimum of F(X), or B is near a C zero of even multiplicity, or neither of these. C C 5 Too many (.GT. 500) function evaluations used. C***REFERENCES L. F. SHAMPINE AND H. A. WATTS, *FZERO, A ROOT-SOLVING C CODE*, SC-TM-70-631, SEPTEMBER 1970. C T. J. DEKKER, *FINDING A ZERO BY MEANS OF SUCCESSIVE C LINEAR INTERPOLATION*, 'CONSTRUCTIVE ASPECTS OF THE C FUNDAMENTAL THEOREM OF ALGEBRA', EDITED BY B. DEJON C P. HENRICI, 1969. C***ROUTINES CALLED R1MACH C***END PROLOGUE FZERO REAL A,ACBS,ACMB,AE,AW,B,C,CMB,ER,FA,FB,FC,FX,FZ,P,Q,R REAL RE,RW,T,TOL,Z INTEGER IC,IFLAG,KOUNT C C ER IS TWO TIMES THE COMPUTER UNIT ROUNDOFF VALUE WHICH IS C DEFINED HERE BY THE FUNCTION R1MACH. C C***FIRST EXECUTABLE STATEMENT FZERO ER = 2.0E0 * R1MACH(4) C C INITIALIZE C Z=R IF(R.LE.AMIN1(B,C).OR.R.GE.AMAX1(B,C)) Z=C RW=AMAX1(RE,ER) AW=AMAX1(AE,0.0) IC=0 T=Z FZ=F(T) FC=FZ T=B FB=F(T) KOUNT=2 IF(SIGN(1.0E0,FZ).EQ.SIGN(1.0E0,FB)) GO TO 1 C=Z GO TO 2 1 IF(Z.EQ.C) GO TO 2 T=C FC=F(T) KOUNT=3 IF(SIGN(1.0E0,FZ).EQ.SIGN(1.0E0,FC)) GO TO 2 B=Z FB=FZ 2 A=C FA=FC ACBS=ABS(B-C) FX=AMAX1(ABS(FB),ABS(FC)) C 3 IF (ABS(FC) .GE. ABS(FB)) GO TO 4 C PERFORM INTERCHANGE A=B FA=FB B=C FB=FC C=A FC=FA C 4 CMB=0.5*(C-B) ACMB=ABS(CMB) TOL=RW*ABS(B)+AW C C TEST STOPPING CRITERION AND FUNCTION COUNT C IF (ACMB .LE. TOL) GO TO 10 IF(FB.EQ.0.E0) GO TO 11 IF(KOUNT.GE.500) GO TO 14 C C CALCULATE NEW ITERATE IMPLICITLY AS B+P/Q C WHERE WE ARRANGE P .GE. 0. C THE IMPLICIT FORM IS USED TO PREVENT OVERFLOW. C P=(B-A)*FB Q=FA-FB IF (P .GE. 0.) GO TO 5 P=-P Q=-Q C C UPDATE A AND CHECK FOR SATISFACTORY REDUCTION C IN THE SIZE OF THE BRACKETING INTERVAL. C IF NOT, PERFORM BISECTION. C 5 A=B FA=FB IC=IC+1 IF (IC .LT. 4) GO TO 6 IF (8.*ACMB .GE. ACBS) GO TO 8 IC=0 ACBS=ACMB C C TEST FOR TOO SMALL A CHANGE C 6 IF (P .GT. ABS(Q)*TOL) GO TO 7 C C INCREMENT BY TOLERANCE C B=B+SIGN(TOL,CMB) GO TO 9 C C ROOT OUGHT TO BE BETWEEN B AND (C+B)/2. C 7 IF (P .GE. CMB*Q) GO TO 8 C C USE SECANT RULE C B=B+P/Q GO TO 9 C C USE BISECTION C 8 B=0.5*(C+B) C C HAVE COMPLETED COMPUTATION FOR NEW ITERATE B C 9 T=B FB=F(T) KOUNT=KOUNT+1 C C DECIDE WHETHER NEXT STEP IS INTERPOLATION OR EXTRAPOLATION C IF (SIGN(1.0,FB) .NE. SIGN(1.0,FC)) GO TO 3 C=A FC=FA GO TO 3 C C C FINISHED. PROCESS RESULTS FOR PROPER SETTING OF IFLAG C 10 IF (SIGN(1.0,FB) .EQ. SIGN(1.0,FC)) GO TO 13 IF (ABS(FB) .GT. FX) GO TO 12 IFLAG = 1 RETURN 11 IFLAG = 2 RETURN 12 IFLAG = 3 RETURN 13 IFLAG = 4 RETURN 14 IFLAG = 5 RETURN END