C #ROOTA V1A 18-MAY-73. C LAST UPDATE: C SUBROUTINE ROOTA (F, A, B, XTOL, FTOL, NTOL, IFLAG, LOG) C C *PURPOSE. C TO PERFORM A MODIFIED REGULA-FALSI SEARCH FOR A ROOT C OF THE USER SUPPLIED FUNCTION "F". C C *PARAMETERS: C C F - THE EXTERNAL NAME OF A USER SUPPLIED FUNCTION FOR C WHICH A ROOT IS REQUIRED. C A - REAL VALUE BEING ONE BOUND ON REGION IN WHICH THE ROOT C IS TO BE FOUND. C B - REAL VALUE BEING THE OTHER BOUND ON THE REGION IN WHICH C THE ROOT IS TO BE FOUND. C XTOL- REAL VALUE BEING THE TOLERABLE ERROR IN THE ROOT POSITION. C FTOL- REAL VALUE BEING THE TOLERABLE DEVIATION OF F(X) FROM C ZERO. C NTOL- INTEGER VALUE BEING MAXIMUM ALLOWED CALLS TO FUNCTION C F(X). C IFLAG- INTEGER VALUE RETURNED GIVING OUTCOME OF SEARCH - SEE C BELOW FOR DETAILS. C LOG - INTEGER VALUE GIVING FORTRAN I/O CHANNEL FOR LOGGING C MESSAGE, OR SET ZERO IF NONE REQUIRED. C THE LOG OUTPUT IS USED FOR ERROR MESSAGES AND SUMMARY OF C THE FINAL RESULT. C C INPUT PARAMETERS: C F, A, B, XTOL, FTOL, NTOL, LOG C OUTPUT PARAMETERS: C A, B, IFLAG C C *METHOD. C THE BOUNDS A AND B ARE REFINED UNTIL EITHER: C (B-A)/2 <= XTOL RETURNS IFLAG=0 C F(X) <= FTOL RETURNS IFLAG=1 C NR OF F(X) CALLS > NTOL RETURNS IFLAG=2 C C IF F(X) HAS SAME SIGN AT ENDPOINTS RETURNS IFLAG=3 AND C IMMEDIATELY RETURNS. C ALTHOUGH THIS IS A RELATIVELY SIMPLE ROOT FINDING C METHOD IT IS USUALLY QUITE EFFICIENT AND HAS THE DECIDED C VIRTUE OF RIGOROUSLY MAINTAINING UPPER AND LOWER BOUNDS FOR C THE ROOT. TYPICALLY THE ROUTINE WILL DO BETTER THAN HALVING THE C RANGE FOR EACH FUNCTION VALUE OBTAINED, SO THIS METHOD C IS TO BE PREFERRED OVER A BINARY SEARCH. C C *ACCURACY. C C *RESTRICTIONS. C THE VALUES OF F(A) AND F(B) MUST BE OF OPPOSITE SIGN. C C *ERROR CONDITIONS. C IFLAG=3 SIGNIFIES THE ABOVE RESTRICTION WAS VIOLATED. C C *NON STANDARD ROUTINES CALLED: C C *TYPICAL TIMES. C C *ORIGIN. C TAKEN FROM CONTE AND DE BOOR: " ELEMENTARY NUMERICAL C ANALYSIS". 1972. BY M.R.MANNING C C *COMMENTS. C C #END. C IFLAG= 0 FA= F(A) SIGNFA= SIGN(1., FA) FB= F(B) CHECK FOR SIGN CHANGE IF (SIGNFA*FB .LE. 0.) GO TO 5 IFLAG= 3 IF (LOG .GT. 0) WRITE(LOG,901) A, B 901 FORMAT(' ROOTA: F(X) HAS SAME SIGN AT ENDPOINTS ',2G14.6) RETURN C 5 W= A FW= FA DO 20 N= 1, NTOL CHECK FOR SUFFICIENTLY SMALL INTERVAL IF (ABS(B-A)/2. .LE. XTOL) GO TO 100 CHECK FOR SUFFICIENTLY SMALL FUNCTION VALUE IF (ABS(FW) .GT. FTOL) GO TO 9 A= W B= W IFLAG= 1 GO TO 100 C 9 W= (FA*B - FB*A) / (FA - FB) PREVFW= SIGN(1., FW) FW= F(W) C CHANGE TO NEW INTERVAL IF (SIGNFA*FW .LT. 0.) GO TO 10 A= W FA= FW IF (FW*PREVFW .GT. 0.) FB= FB/2. GO TO 20 C 10 B= W FB= FW IF (FW*PREVFW .GT. 0.) FA= FA/2. 20 CONTINUE IFLAG= 2 IF (LOG .GT. 0) WRITE(LOG, 902) NTOL 902 FORMAT(' ROOTA: NO CONVERGENCE IN ',I5,' ITERATIONS') C 100 IF (LOG .GT. 0) WRITE(LOG, 903) A, B 903 FORMAT(' ROOTA: ROOT IN RANGE ',G14.6,' TO ', G14.6) RETURN C END