C #ROOTB V1A 28-JUN-73. C LAST UPDATE: C SUBROUTINE ROOTB(FUN, A, B, XTOL, FTOL, NTOL, IFLAG, LOG) C C *PURPOSE. C TO FIND A ROOT OF THE USER SUPPLIED FUNCTION "FUN(X)" IN C THE RANGE "A" TO "B". UNLIKE "ROOTA" THIS ROUTINE DOES NOT C REQUIRE THAT "FUN" HAVE OPPOSITE SIGNS AT THE ENDPOINTS, AND C IF THE SAME SIGN IS FOUND "ROOTB" WILL SET UP A SEARCH WHICH C EITHER FINDS A SIGN CHANGE AND THEN A ROOT OR ELSE A STATIONARY C POINT. C C *PARAMETERS: C C FUN - THE EXTERNAL NAME OF THE USER'S FUNCTION FOR WHICH A C ROOT IS REQUIRED. C A - REAL VALUE OF ONE ENDPOINT OF RANGE TO BE USED. C B - REAL VALUE OF OTHER ENDPOINT OF RANGE TO BE USED. C XTOL - REAL VALUE OF TOLERABLE ERROR IN POSITION OF ROOT. C FTOL - REAL VALUE OF TOLERABLE ERROR IN FUNCTION VALUE AT AN C ESTIMATED ROOT POSITION. C NTOL - INTEGER MAXIMUM NUMBER OF EVALUATIONS OF "FUN(X)" TO BE C CARRIED OUT. C IFLAG- INTEGER VARIABLE RETURNING THE OUTCOME OF THE SEARCH AS: C 0 - INTERVAL (A,B) REDUCED TO LESS THAN 2*XTOL AFTER C SIGN CHANGE FOUND. I.E. ROOT FOUND! C 1 - FUNCTION VALUE LESS THAN "FTOL" FOUND AT X=A=B AFTER C SIGN CHANGE FOUND. I.E. ROOT FOUND! C 2 - MAX NR OF "FUN" EVALUATIONS DONE AND NO CONVERGENCE C BUT A SIGN CHANGE HAS BEEN FOUND. C 3 - INTERVAL (A,B) REDUCED TO LESS THAN 2*XTOL BEFORE C ANY SIGN CHANGE FOUND. I.E. LOCAL STATIONARY POINT. C 4 - MAX NR OF "FUN" EVALUATIONS DONE AND NO CONVERGENCE C AND ALSO NO SIGN CHANGE FOUND!! C LOG - INTEGER VALUE EITHER ZERO MEANING NO LOGGING MESSAGES C TO BE WRITTEN OUT, OR ELSE IS THE FORTRAN I/O CHANNEL C TO WRITE THESE MESSAGES TO. C C NOTES. THE INTERVAL ENDPOINTS "A" AND "B" ARE BROUGHT CLOSER C TOGETHER AND SHOULD ALWAYS BRACKET THE TARGET POINT. C THE LOGGING MESSAGES ARE ALL CLEARLY IDENTIFIED AS COMING C FROM "ROOTB" AND GIVE INDICATION OF NO CONVERGENCE, NO SIGN C CHANGE FOUND IF APPROPRIATE AS WELL AS SUMMARISING THE INITIAL C AND FINAL INTERVALS AND THE "FUN" VALUES AT THE ENDS. C C INPUT PARAMETERS: C FUN, A, B, XTOL, FTOL, NTOL, LOG C OUTPUT PARAMETERS: C A, B, IFLAG C C *METHOD. C THE ROOT SEARCHING METHOD IS THE MODIFIED REGULA-FALSI C METHOD AS USED IN "ROOTA", SEE FOR EXAMPLE CONTE AND DE BOOR, C "ELEMENTARY NUMERICAL ANALYSIS". THIS METHOD HAS THE NICE C PROPERTIES OF ENSURING THE ROOT POSITION IS BRACKETED AND C GENERALLY CONVERGING FASTER THAN A BINARY BISECTION C METHOD. HOWEVER THIS SECTION OF THE ROUTINE CAN ONLY BE USED C WHEN A SIGN CHANGE HAS BEEN FOUND. C IF THE VALUES OF "FUN" AT THE ENDS OF THE INITIAL INTERVAL C ARE OF THE SAME SIGN "ROOTB" SETS UP A "GOLDEN-SECTION" SEARCH C FOR AN EXTREMUM POINT. THUS IF "FUN" IS FOUND POSITIVE THE C SEARCH IS FOR A MINIMUM AND IF "FUN" IS NEGATIVE FOR A MAXIMUM. C THE GOLDEN-SECTION SEARCH ( SEE FOR EXAMPLE J.KOWALIK AND C M.OSBORNE, "METHODS FOR UNCONSTRAINED OPTIMIZATION PROBLEMS") C HAS THE NICE PROPERTY OF MAINTAINING A BRACKETING INTERVAL C ON THE EXTREMUM. THIS SEARCH WILL BE STOPPED AS SOON AS A SIGN C CHANGE IS FOUND AND THE ROOT-SEARCH STARTED UP, HOWEVER C IF NO SIGN CHANGE CAN BE FOUND THIS SEARCH WILL PRODUCE A VALUE C FOR THE STATIONARY POINT. C IF THERE ARE TWO ROOTS IN THE INTERVAL THE ONE CLOSEST C TO "A" WILL BE FOUND. THE SIMPLEST TYPE OF FUNCTION WHICH WILL C POSSIBLY CAUSE TROUBLE IS ONE WITH TWO ROOTS AND TWO EXTREMA C IN THE RANGE "A" TO "B". IN THIS CASE THE ROUTINE MAY FIND C ONE OF THE EXTREMA INSTEAD OF ONE OF THE ROOTS. C NOTE THAT THE ROUTINE DOES NOT STOP WHEN ENCOUNTERING C THE CASE OF "FUN" LESS THAN "FTOL" BEFORE A SIGN CHANGE FOUND. C THIS IS TO PREVENT ASYMPTOTIC BEHAVIOUR BEING MISTAKEN FOR C AN INDICATION OF A ROOT. C ALTHOUGH THE METHOD CAN BE USED TO FIND MAXIMA OR MINIMA C IN A GIVEN RANGE IT WILL GENERALLY NOT BE VERY EFFICIENT IN C THIS RESPECT UNLESS THE FUNCTION IS SHARPLY VARYING. C C *ACCURACY. C C *RESTRICTIONS. C SEE ABOVE FOR PROBLEM CASES. C C *ERROR CONDITIONS. C SEE "IFLAG" ABOVE. C C *NON STANDARD ROUTINES CALLED. C "FUN" - USER'S ROUTINE. C C *TYPICAL TIMES. C C *ORIGIN. C M.R.MANNING WITH HELP FROM BOTH OF THE ABOVE REFERENCES. C C *COMMENTS. C C #END. C C NOTE THIS IS A HORRIBLY UNSTRUCTURED PROGRAM - NASTY GO-TO'S C EVERYWHERE - BUT IT SEEMS THE ONLY EFFICIENT WAY OF DOING IT C AND THE COMMENTS SHOULD EXPLAIN THE MESS. C INTEGER CYCLE C C************************************************************ C C INITIALIZATION SECTION: C SETS UP INITIAL VALUES AND CHECKS SIGN CHANGE OVER INITIAL C RANGE. C IFLAG= 3 INOV= 4 FA= FUN(A) SIGNFA= SIGN(1., FA) FB= FUN(B) C REMEMBER INTIAL VALUES FOR LATER LOG OUTPUT. SA= A SB= B SFA= FA SFB= FB C NEVAL= 2 C CHECK CHANGE OF SIGN IF ( SIGNFA .NE. SIGN(1., FB) ) GO TO 200 C C************************************************************ C C SET UP GOLDEN-SECTION SEARCH FOR EXTREMUM IN THE INTERVAL C A TO B. I.E. LOOKING FOR MINIMUM IF FUN +VE, AND MAXIMUM C IF FUN IS -VE. C GFAC= (3. - SQRT(5.) )/2. C TEMPORARY SWITCH OF "CYCLE" LABEL. ASSIGN 50 TO CYCLE GO TO 160 C NOW RESET "CYCLE" TO GOLDEN-SECTION SEARCH. 50 ASSIGN 100 TO CYCLE GO TO 180 C THAT FINISHES INITIALIZATION - NOW EVERYTHING RUNS. C C************************************************************ C C CYCLE SECTION FOR GOLDEN-SECTION SEARCH. SHIFTS ONE C END-POINT AND RECALCULATES ONE NEW FUNCTION VALUE. IF A C SIGN CHANGE IS FOUND EXITS TO ROOT SEARCH SECTION BELOW C AT STATEMENT 200. C 100 IF (SIGNFA .NE. SIGN(1., F2-F1) ) GO TO 120 C NOTE THAT THAT TEST IS LOOKING FOR A MAX OR MINIMA DEPENDING C ON THE SIGN OF "FA". C C SHIFT DOWN SECTION. C B= X2 FB= F2 X2= X1 F2= F1 C C EVALUATE F1 AT X1 , CHECK SIGN, AND CYCLE. C 160 X1= A + GFAC*(B - A) F1= FUN(X1) IF (SIGNFA .EQ. SIGN(1., F1) ) GO TO 320 C C EUREKA! - FOUND A SIGN CHANGE START UP ROOT-SEARCH. C B= X1 FB= F1 GO TO 200 C C SHIFT UP SECTION. C 120 A= X1 FA= F1 X1= X2 F1= F2 C C EVALUATE F2 AT X2 , CHECK SIGN, AND CYCLE. C 180 X2= B - GFAC*(B - A) F2= FUN(X2) IF (SIGNFA .EQ. SIGN(1., F2) ) GO TO 320 C C EUREKA! - FOUND SIGN CHANGE START UP ROOT-SEARCH. C B= X2 FB= F2 C GO TO 200 C C************************************************************* C C START UP THE MODIFIED REGULA-FALSI SEARCH FOR A ROOT IN C THE RANGE "A" TO "B" NOW THAT A SIGN CHANGE HAS BEEN C ESTABLISHED OVER THAT INTERVAL. C 200 W= A FW= FA SIGNFW= SIGNFA IFLAG= 0 INOV= 2 C SET "CYCLE" TO ROOT-SEARCH SECTION. ASSIGN 250 TO CYCLE GO TO 300 C C************************************************************ C C CYCLE TO HERE IN ROOT-SEARCH PHASE. MAKE A NEW ESTIMATE C OF THE ROOT POSITION, MOVE TO THAT POINT AND SHIFT ONE C ENDPOINT. ALSO FUDGE THE CHORD TO IMPROVE THE REGULA FALSI C METHOD. C 250 W= (FA*B - FB*A)/(FA - FB) PREVFW= SIGNFW FW= FUN(W) SIGNFW= SIGN(1., FW) IF (SIGNFA .NE. SIGNFW) GO TO 270 C NEW POINT HAS SAME SIGN AS ENDPOINT A - SO MOVE A. A= W FA= FW IF ( SIGNFW .EQ. PREVFW) FB=FB/2. GO TO 300 C C NEW POINT HAS SAME SIGN AS ENDPOINT B - SO MOVE B. 270 B= W FB= FW IF (SIGNFW .EQ. PREVFW) FA= FA/2. C GO TO 300 C C************************************************************* C C THIS SECTION CHECKS ALL THE CONDITIONS FOR TERMINATING C THE ABOVE PROCESSES. IT IS ENTERED AFTER EVERY NEW C FUNCTION EVALUATION (EXCEPT THE FIRST 2) AND WILL RETURN C TO THE LABEL "CYCLE" SOMEWHERE ABOVE OR ELSE EXIT FROM C THE SUBROUTINE AFTER CLEAN UP AND OPTIONAL LOG MESSAGES. C 300 IF (ABS(FW) .GT. FTOL) GO TO 320 IFLAG= 1 A= W B= W GO TO 400 C C BEFORE A SIGN CHANGE IS DETECTED THE ABOVE TEST IS NOT C APPLIED, SO THE CHECK ENTRY POINT IS HERE. C 320 IF (ABS(B - A)/2. .LT. XTOL) GO TO 400 NEVAL= NEVAL + 1 IF (NEVAL .LT. NTOL) GO TO CYCLE C OOPS! - TOO MANY FUNCTION VALUES, INDICATE LACK OF CONVERGENCE IFLAG= INOV IF (LOG .GT. 0) WRITE(LOG, 902) NTOL 902 FORMAT(' ROOTB:***NO CONVERGENCE AFTER ',I6,' EVALUATIONS***') C C GENERAL EXIT POINT: C CHECK IF SIGN CHANGE FOUND AND WARN IF NOT, C PRINT STATISTICS IF WANTED. C 400 IF ( INOV .EQ. 4 .AND. LOG .GT. 0) WRITE(LOG, 901) 901 FORMAT(' ROOTB:***WARNING NO SIGN CHANGE FOUND***') IF ( LOG .LE. 0) RETURN WRITE(LOG, 903) 903 FORMAT(' ROOTB: RANGE REDUCED FROM:') WRITE(LOG, 900) SA, SFA, SB, SFB 900 FORMAT(' ROOTB: A ',G13.6,' F(A) ',G13.6,' B ',G13.6 1 ,' F(B) ',G13.6) WRITE(LOG, 904) 904 FORMAT(' ROOTB: TO:') WRITE(LOG, 900) A, FA, B, FB RETURN C END