*DECK DDZRO SUBROUTINE DDZRO (AE, F, H, N, NQ, IROOT, RE, T, YH, UROUND, B, C, 8 FB, FC, Y) C***BEGIN PROLOGUE DDZRO C***SUBSIDIARY C***PURPOSE DDZRO SEARCHES FOR A ZERO OF A FUNCTION F(N, T, Y, IROOT) C BETWEEN THE GIVEN VALUES B AND C UNTIL THE WIDTH OF THE C INTERVAL (B, C) HAS COLLAPSED TO WITHIN A TOLERANCE C SPECIFIED BY THE STOPPING CRITERION, C ABS(B - C) .LE. 2.*(RW*ABS(B) + AE). C***LIBRARY SLATEC (SDRIVE) C***TYPE DOUBLE PRECISION (SDZRO-S, DDZRO-D, CDZRO-C) C***AUTHOR KAHANER, D. K., (NIST) C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C SUTHERLAND, C. D., (LANL) C MAIL STOP D466 C LOS ALAMOS NATIONAL LABORATORY C LOS ALAMOS, NM 87545 C***DESCRIPTION C C THIS IS A SPECIAL PURPOSE VERSION OF ZEROIN, MODIFIED FOR USE WITH C THE DDRIV PACKAGE. C C SANDIA MATHEMATICAL PROGRAM LIBRARY C MATHEMATICAL COMPUTING SERVICES DIVISION 5422 C SANDIA LABORATORIES C P. O. BOX 5800 C ALBUQUERQUE, NEW MEXICO 87115 C CONTROL DATA 6600 VERSION 4.5, 1 NOVEMBER 1971 C C PARAMETERS C F - NAME OF THE EXTERNAL FUNCTION, WHICH RETURNS A C DOUBLE PRECISION RESULT. THIS NAME MUST BE IN AN C EXTERNAL STATEMENT IN THE CALLING PROGRAM. 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 - THE OTHER END OF THE INTERVAL (B, 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 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***REFERENCES L. F. SHAMPINE AND H. A. WATTS, ZEROIN, A ROOT-SOLVING C ROUTINE, SC-TM-70-631, SEPT 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 AND P. HENRICI, 1969. C***ROUTINES CALLED DDNTP C***REVISION HISTORY (YYMMDD) C 790601 DATE WRITTEN C 900329 INITIAL SUBMISSION TO SLATEC. C***END PROLOGUE DDZRO INTEGER IC, IROOT, KOUNT, N, NQ DOUBLE PRECISION A, ACBS, ACMB, AE, B, C, CMB, ER, F, FA, FB, FC, 8 H, P, Q, RE, RW, T, TOL, UROUND, Y(*), YH(N,*) C***FIRST EXECUTABLE STATEMENT DDZRO ER = 4.D0*UROUND RW = MAX(RE, ER) IC = 0 ACBS = ABS(B - C) A = C FA = FC KOUNT = 0 C PERFORM INTERCHANGE 10 IF (ABS(FC) .LT. ABS(FB)) THEN A = B FA = FB B = C FB = FC C = A FC = FA END IF CMB = 0.5D0*(C - B) ACMB = ABS(CMB) TOL = RW*ABS(B) + AE C TEST STOPPING CRITERION IF (ACMB .LE. TOL) RETURN IF (KOUNT .GT. 50) RETURN C CALCULATE NEW ITERATE IMPLICITLY AS C B + P/Q, WHERE WE ARRANGE P .GE. 0. C THE IMPLICIT FORM IS USED TO PREVENT OVERFLOW. P = (B - A)*FB Q = FA - FB IF (P .LT. 0.D0) THEN P = -P Q = -Q END IF C UPDATE A AND CHECK FOR SATISFACTORY REDUCTION C IN THE SIZE OF OUR BOUNDING INTERVAL. A = B FA = FB IC = IC + 1 IF (IC .GE. 4) THEN IF (8.D0*ACMB .GE. ACBS) THEN C BISECT B = 0.5D0*(C + B) GO TO 20 END IF IC = 0 END IF ACBS = ACMB C TEST FOR TOO SMALL A CHANGE IF (P .LE. ABS(Q)*TOL) THEN C INCREMENT BY TOLERANCE B = B + SIGN(TOL, CMB) C ROOT OUGHT TO BE BETWEEN C B AND (C + B)/2. ELSE IF (P .LT. CMB*Q) THEN C INTERPOLATE B = B + P/Q ELSE C BISECT B = 0.5D0*(C + B) END IF C HAVE COMPLETED COMPUTATION C FOR NEW ITERATE B. 20 CALL DDNTP (H, 0, N, NQ, T, B, YH, Y) FB = F(N, B, Y, IROOT) IF (N .EQ. 0) RETURN IF (FB .EQ. 0.D0) RETURN KOUNT = KOUNT + 1 C C DECIDE WHETHER NEXT STEP IS INTERPOLATION OR EXTRAPOLATION C IF (SIGN(1.0D0, FB) .EQ. SIGN(1.0D0, FC)) THEN C = A FC = FA END IF GO TO 10 END *DECK DDNTP SUBROUTINE DDNTP (H, K, N, NQ, T, TOUT, YH, Y) C***BEGIN PROLOGUE DDNTP C***SUBSIDIARY C***PURPOSE SUBROUTINE DDNTP INTERPOLATES THE K-TH DERIVATIVE OF Y AT C TOUT, USING THE DATA IN THE YH ARRAY. IF K HAS A VALUE C GREATER THAN NQ, THE NQ-TH DERIVATIVE IS CALCULATED. C***LIBRARY SLATEC (SDRIVE) C***TYPE DOUBLE PRECISION (SDNTP-S, DDNTP-D, CDNTP-C) C***AUTHOR KAHANER, D. K., (NIST) C NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY C GAITHERSBURG, MD 20899 C SUTHERLAND, C. D., (LANL) C MAIL STOP D466 C LOS ALAMOS NATIONAL LABORATORY C LOS ALAMOS, NM 87545 C***ROUTINES CALLED (NONE) C***REVISION HISTORY (YYMMDD) C 790601 DATE WRITTEN C 900329 INITIAL SUBMISSION TO SLATEC. C***END PROLOGUE DDNTP INTEGER I, J, JJ, K, KK, KUSED, N, NQ DOUBLE PRECISION FACTOR, H, R, T, TOUT, Y(*), YH(N,*) C***FIRST EXECUTABLE STATEMENT DDNTP IF (K .EQ. 0) THEN DO 10 I = 1,N 10 Y(I) = YH(I,NQ+1) R = ((TOUT - T)/H) DO 20 JJ = 1,NQ J = NQ + 1 - JJ DO 20 I = 1,N 20 Y(I) = YH(I,J) + R*Y(I) ELSE KUSED = MIN(K, NQ) FACTOR = 1.D0 DO 40 KK = 1,KUSED 40 FACTOR = FACTOR*(NQ+1-KK) DO 50 I = 1,N 50 Y(I) = FACTOR*YH(I,NQ+1) R = ((TOUT - T)/H) DO 80 JJ = KUSED+1,NQ J = KUSED + 1 + NQ - JJ FACTOR = 1.D0 DO 60 KK = 1,KUSED 60 FACTOR = FACTOR*(J-KK) DO 70 I = 1,N 70 Y(I) = FACTOR*YH(I,J) + R*Y(I) 80 CONTINUE DO 100 I = 1,N 100 Y(I) = Y(I)*H**(-KUSED) END IF RETURN END