DIMENSION Y(101) C C C MAIN PROGRAM AND USER WRITTEN SUBROUTINES FOR THE FUNCTION C F WHOSE KTH COMPONENT IS C EXP(COS(K*(X(1)+X(2)+...+X(N)))) . C C N=3 ARCTOL=1.0E-4 EPS=1.0E-8 CALL SECOND(TIME1) IFLAG=0 CALL FIXPT(N,Y,ARCTOL,EPS,ARCLEN,NFE,IFLAG) NP1=N+1 CALL SECOND(TIME2) WRITE (6,30) ARCLEN,IFLAG,(Y(L),L=1,NP1) 30 FORMAT(/11H ARC LENGTH,1PE18.9,4X,4HCODE,I3/(1X,7E17.9)) WRITE (6,40) NFE 40 FORMAT(/I8,21H JACOBIAN EVALUATIONS/) EXTIME=TIME2-TIME1 WRITE (6,50) EXTIME 50 FORMAT(22H EXECUTION TIME(SEC) =,F8.2/1H1) 60 CONTINUE STOP END SUBROUTINE SECOND(T) C C **** THIS IS A DUMMY ROUTINE **** C C T SHOULD BE SET EQUAL TO A COUNTER COUNTING SECONDS C T = 0.0 RETURN END SUBROUTINE F(X,V) C C C SUBROUTINE TO EVALUATE THE FUNCTION WHOSE JTH COMPONENT IS C EXP(COS(J*(X(1)+...+X(N)))) . C C ON RETURN V CONTAINS F(X) . C C REAL X(100),V(100) COMMON /FIXEDP/ YPOLD(101),A(100),NFE,N,NP1,IFLAG SUMX=0.0 DO 20 J=1,N 20 SUMX=SUMX+X(J) DO 30 K=1,N 30 V(K)=EXP(COS(FLOAT(K)*SUMX)) RETURN END SUBROUTINE FJAC(X,V,K) C C C EVALUATES THE JACOBIAN MATRIX OF THE FUNCTION WHOSE JTH C COMPONENT IS C EXP(COS(J*(X(1)+...+X(N)))) . C C ON RETURN V CONTAINS THE KTH COLUMN OF THE JACOBIAN MATRIX C DF(X) . C C C REAL X(100),V(100) COMMON /FIXEDP/ YPOLD(101),A(100),NFE,N,NP1,IFLAG IF (K .GT. 1) RETURN SUMX=0.0 DO 30 J=1,N 30 SUMX=SUMX+X(J) 40 DO 60 J=1,N FJ=FLOAT(J) ECKX=EXP(COS(FJ*SUMX)) 60 V(J)=-FJ*SIN(FJ*SUMX)*ECKX RETURN C C C END OF TEST PROGRAMS. C C **************************************************************** C END SUBROUTINE FIXPT(N,Y,ARCTOL,EPS,ARCLEN,NFE,IFLAG) C C SUBROUTINE FIXPT FINDS A FIXED POINT OR ZERO OF THE C N-DIMENSIONAL VECTOR FUNCTION F(X). FOR THE FIXED POINT C PROBLEM F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL INTO C ITSELF. THE EQUATION X = F(X) IS SOLVED BY C FOLLOWING THE ZERO CURVE OF THE HOMOTOPY MAP C C LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A) , C C STARTING FROM LAMBDA = 0, X = A. THE CURVE IS PARAMETERIZED C BY ARC LENGTH S, AND IS FOLLOWED BY SOLVING THE ORDINARY C DIFFERENTIAL EQUATION D(HOMOTOPY MAP)/DS = 0 FOR C Y(S) = (LAMBDA(S), X(S)). C C FOR THE ZERO FINDING PROBLEM F(X) IS ASSUMED TO BE A C2 MAP C SUCH THAT FOR SOME R .GT. 0, X*F(X) .GE. 0 WHENEVER NORM(X) = R. C THE EQUATION F(X) = 0 IS SOLVED BY FOLLOWING THE ZERO CURVE C OF THE HOMOTOPY MAP C C LAMBDA*F(X) + (1 - LAMBDA)*(X - A) C C EMANATING FROM LAMBDA = 0, X = A. C C A MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS. C C THIS CODE IS BASED ON THE ALGORITHM IN L. T. WATSON, A C GLOBALLY CONVERGENT ALGORITHM FOR COMPUTING FIXED POINTS OF C C2 MAPS, APPL. MATH. COMPUT., 1978. C C C THE USER MUST SUPPLY A SUBROUTINE F(X,V) WHICH EVALUATES C F(X) AT X AND RETURNS THE VECTOR F(X) IN V, AND A SUBROUTINE C FJAC(X,V,K) WHICH RETURNS IN V THE KTH COLUMN OF THE C JACOBIAN MATRIX OF F(X) EVALUATED AT X. FIXPT DIRECTLY OR C INDIRECTLY USES THE SUBROUTINES STEP , INTRP , ROOT , FODE , C F , FJAC , DCPOSE , AND THE FUNCTION INNER . FIXPT AND C FODE CONTAIN DIMENSION STATEMENTS WHICH LIMIT N TO 100. C STEP AND ROOT CONTAIN MACHINE DEPENDENT CONSTANTS. NO C OTHER MODIFICATIONS BY THE USER ARE REQUIRED. C C C ON INPUT: C C N IS THE DIMENSION OF X AND F(X). C C Y IS AN ARRRAY OF LENGTH N + 1. (Y(2),...,Y(N+1)) = A IS THE C STARTING POINT FOR THE ZERO CURVE. C C ARCTOL IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN C FOLLOWING THE ZERO CURVE. IF ARCTOL .LE. 0.0 ON INPUT C IT IS RESET TO .5*SQRT(EPS). NORMALLY ARCTOL SHOULD C BE CONSIDERABLY LARGER THAN EPS. C C EPS IS THE LOCAL ERROR ALLOWED THE ODE SOLVER WHEN VERY C NEAR THE FIXED POINT(ZERO). EPS IS APPROXIMATELY THE C ABSOLUTE ERROR IN THE COMPUTED FIXED POINT(ZERO). C C IFLAG CAN BE -1, 0, 2, OR 3. IFLAG SHOULD BE 0 ON THE FIRST C CALL TO FIXPT FOR THE PROBLEM X=F(X), AND -1 FOR THE C PROBLEM F(X)=0. IN CERTAIN SITUATIONS IFLAG IS SET TO C 2 OR 3 BY FIXPT, AND FIXPT CAN BE CALLED AGAIN WITHOUT C CHANGING IFLAG. C C Y, ARCTOL, EPS, ARCLEN, NFE, AND IFLAG SHOULD ALL BE C VARIABLES IN THE CALLING PROGRAM. C C C ON OUTPUT: C C N IS UNCHANGED. C C Y(1) = LAMBDA, (Y(2),...,Y(N+1)) = X, AND Y IS AN APPROXIMATE C ZERO OF THE HOMOTOPY MAP. NORMALLY LAMBDA = 1 AND X IS A C FIXED POINT(ZERO) OF F(X). IN ABNORMAL SITUATIONS LAMBDA C MAY ONLY BE NEAR 1 AND X IS NEAR A FIXED POINT(ZERO). C C ARCTOL = EPS AFTER A NORMAL RETURN (IFLAG = 1). C C EPS IS UNCHANGED AFTER A NORMAL RETURN (IFLAG = 1). IT IS C INCREASED TO AN APPROPRIATE VALUE ON THE RETURN IFLAG = 2. C C ARCLEN IS THE LENGTH OF THE PATH FOLLOWED. C C NFE IS THE NUMBER OF FUNCTION EVALUATIONS (= NUMBER OF C JACOBIAN EVALUATIONS). C C IFLAG = C -1 CAUSES FIXPT TO INITIALIZE EVERYTHING FOR THE PROBLEM C F(X) = 0 (USE ON FIRST CALL). C C 0 CAUSES FIXPT TO INITIALIZE EVERYTHING FOR THE PROBLEM C X = F(X) (USE ON FIRST CALL). C C 1 NORMAL RETURN. C C 2 SPECIFIED ERROR TOLERANCE CANNOT BE MET. EPS HAS BEEN C INCREASED TO A SUITABLE VALUE. TO CONTINUE, JUST CALL C FIXPT AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 3 STEP HAS BEEN CALLED 1000 TIMES. TO CONTINUE, CALL C FIXPT AGAIN WITHOUT CHANGING ANY PARAMETERS. C C 4 JACOBIAN MATRIX DOES NOT HAVE FULL RANK. THE ALGORITHM C HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE C FOLLOWED ANY FURTHER). C C 5 EPS (OR ARCTOL ) IS TOO LARGE. THE PROBLEM SHOULD BE C RESTARTED BY CALLING FIXPT WITH A SMALLER EPS (OR C ARCTOL ) AND IFLAG = 0 (-1). C C 6 I - DF(X) IS NEARLY SINGULAR AT THE FIXED POINT(DF(X) IS C NEARLY SINGULAR AT THE ZERO). ANSWER MAY NOT BE C ACCURATE. C C 7 ILLEGAL INPUT PARAMETERS, A FATAL ERROR. C C ***** ARRAY DECLARATIONS. ***** C C THE FOLLOWING ARRAYS ARE FOR A MAXIMUM OF N = 100 VARIABLES. C TO HANDLE UP TO N VARIABLES, CHANGE EACH 100 AND 101 TO N C AND N + 1 RESPECTIVELY. C C ARRAYS NEEDED BY THE ODE SUBROUTINE STEP . REAL Y(101),WT(101),PHI(101,16),P(101),YP(101),PSI(12) C COMMON /FIXEDP/ YPOLD(101),A(100),NFEC,NC,NP1,IFLAGC C C ***** END OF DIMENSIONAL INFORMATION. ***** C EXTERNAL FODE LOGICAL START,CRASH,ST99 C C LIMITD IS AN UPPER BOUND ON THE NUMBER OF STEPS. IT MAY BE C CHANGED BY CHANGING THE FOLLOWING DATA STATEMENT: DATA LIMITD/1000/ C C C : : : : : : : : : : : : : : : : : : : : : IF (N .LE. 0 .OR. EPS .LE. 0.0) IFLAG=7 IF (IFLAG .LE. 0) GO TO 10 IF (IFLAG .EQ. 2) GO TO 35 IF (IFLAG .EQ. 3) GO TO 30 C ONLY VALID INPUT FOR IFLAG IS -1, 0, 2, 3. IFLAG=7 RETURN C C ***** INITIALIZATION BLOCK. ***** C 10 ARCLEN=0.0 S=0.0 IF (ARCTOL .LE. 0.0) ARCTOL=.5*SQRT(EPS) NFEC=0 NC=N IFLAGC=IFLAG NP1=N+1 SQNP1=SQRT(FLOAT(NP1)) CURTOL=.01/SQNP1 ST99=.FALSE. START=.TRUE. CRASH=.FALSE. H=.1 EPSSTP=ARCTOL C SET INITIAL CONDITIONS FOR ORDINARY DIFFERENTIAL EQUATION. YPOLD(1)=1.0 Y(1)=0.0 WT(1)=1.0 DO 20 J=2,NP1 YPOLD(J)=0.0 A(J-1)=Y(J) WT(J)=1.0 20 CONTINUE 30 LIMIT=LIMITD C C ***** END OF INITIALIZATION BLOCK. ***** C C C ***** MAIN LOOP. ***** C 35 DO 150 ITER=1,LIMIT IF (Y(1) .GE. 0.0) GO TO 50 40 IFLAG=5 RETURN 50 IF (S .LE. 7.0*SQNP1) GO TO 80 C ARC LENGTH IS GETTING TOO LONG, THE PROBLEM WILL BE C RESTARTED WITH A DIFFERENT A VECTOR. ARCLEN=ARCLEN+S S=0.0 60 START=.TRUE. CRASH=.FALSE. C COMPUTE A NEW A VECTOR. CALL F(Y(2),YP) DO 70 JW=1,N AOLD=A(JW) IF (IFLAGC .LT. 0) + A(JW)=Y(1)*YP(JW)/(1.0 - Y(1)) + Y(JW+1) IF (IFLAGC .EQ. 0) + A(JW)=(Y(JW+1) - Y(1)*YP(JW))/(1.0 - Y(1)) IF (ABS(A(JW)-AOLD) .GT. .95) GO TO 40 70 CONTINUE GO TO 100 80 IF (Y(1) .LE. .99 .OR. ST99) GO TO 100 C WHEN LAMBDA REACHES .99, THE PROBLEM WILL BE RESTARTED WITH C A NEW A VECTOR. 90 ST99=.TRUE. EPSSTP=EPS ARCTOL=EPS GO TO 60 C C TAKE A STEP ALONG THE CURVE. 100 CALL STEP(S,Y,FODE,NP1,H,EPSSTP,WT,START,HOLD,K,KOLD,CRASH, + PHI,P,YP,PSI) NFE=NFEC C CHECK IF THE STEP WAS SUCCESSFUL. IF (IFLAGC .NE. 4) GO TO 120 IFLAG=4 RETURN 120 IF (.NOT. CRASH) GO TO 130 C RETURN CODE FOR ERROR TOLERANCE TOO SMALL. IFLAG=2 C CHANGE ERROR TOLERANCES. EPS=EPSSTP IF (ARCTOL .LT. EPS) ARCTOL=EPS C CHANGE LIMIT ON NUMBER OF ITERATIONS. LIMS=LIMIT-ITER GO TO 220 C 130 EPSSTP=ARCTOL IF (ABS(YP(1)) .LE. CURTOL) EPSSTP=EPS IF (Y(1) .LT. 1.0) GO TO 150 IF (ST99) GO TO 160 C C IF LAMBDA .GE. 1.0 BUT THE PROBLEM HAS NOT BEEN RESTARTED C WITH A NEW A VECTOR, BACK UP AND RESTART. C S99=S-.5*HOLD C GET AN APPROXIMATE ZERO Y(S) WITH Y(1)=LAMBDA .LT. 1.0 . 135 CALL INTRP(S,Y,S99,WT,P,NP1,KOLD,PHI,PSI) IF (WT(1) .LT. 1.0) GO TO 140 S99=.5*(S-HOLD+S99) GO TO 135 C 140 DO 144 JUDY=1,NP1 Y(JUDY)=WT(JUDY) YPOLD(JUDY)=P(JUDY) 144 WT(JUDY)=1.0 S=S99 GO TO 90 C 150 CONTINUE C C ***** END OF MAIN LOOP. ***** C C LAMBDA HAS NOT REACHED 1 IN 1000 STEPS. IFLAG=3 RETURN C C C USE INVERSE INTERPOLATION TO GET THE ANSWER AT LAMBDA = 1.0 . C 160 SA=S-HOLD SB=S LCODE=1 170 CALL ROOT(SOUT,Y1SOUT,SA,SB,EPS,EPS,LCODE) C ROOT FINDS S SUCH THAT Y(1)(S) = LAMBDA = 1 . IF (LCODE .GT. 0) GO TO 190 CALL INTRP(S,Y,SOUT,WT,P,NP1,KOLD,PHI,PSI) Y1SOUT=WT(1)-1.0 GO TO 170 190 IFLAG=1 C SET IFLAG = 6 IF ROOT COULD NOT GET LAMBDA = 1.0 . IF (LCODE .GT. 2) IFLAG=6 ARCLEN=ARCLEN+SA C LAMBDA(SA) = 1.0 . CALL INTRP(S,Y,SA,WT,P,NP1,KOLD,PHI,PSI) C DO 210 J=1,NP1 210 Y(J)=WT(J) RETURN 220 LIMIT=LIMS RETURN END SUBROUTINE FODE(S,Y,YP) C C SUBROUTINE FODE IS USED BY SUBROUTINE STEP TO SPECIFY THE C ORDINARY DIFFERENTIAL EQUATION DY/DS = G(S,Y) , WHOSE SOLUTION C IS THE ZERO CURVE OF THE HOMOTOPY MAP. S = ARC LENGTH, C YP = DY/DS, AND Y(S) = (LAMBDA(S), X(S)) . C C ***** ARRAY DECLARATIONS. ***** C C THE FOLLOWING ARRAYS ARE FOR A MAXIMUM OF N = 100 VARIABLES. C TO HANDLE UP TO N VARIABLES, CHANGE EACH 100 AND 101 TO N C AND N + 1 RESPECTIVELY. C REAL Y(1),YP(1) C ARRAYS FOR COMPUTING THE JACOBIAN MATRIX AND ITS KERNAL. REAL QR(100,101),ALPHA(100),TZ(101) INTEGER PIVOT(101) C REAL INNER COMMON /FIXEDP/ YPOLD(101),A(100),NFE,N,NP1,IFLAG NDIM=100 C C ***** END OF DIMENSIONAL INFORMATION. ***** C NFE=NFE+1 C NFE CONTAINS THE NUMBER OF JACOBIAN EVALUATIONS. C * * * * * * * * * * * * * * * * * C C COMPUTE THE JACOBIAN MATRIX, STORE IT IN QR. C C QR = ( A - F(X), I - LAMBDA*DF(X) ) . C CALL F(Y(2),TZ) IF (IFLAG .LT. 0) GO TO 140 DO 100 J=1,N 100 QR(J,1)=A(J)-TZ(J) DO 120 K=1,N CALL FJAC(Y(2),TZ,K) KP1=K+1 DO 110 J=1,N 110 QR(J,KP1)=-Y(1)*TZ(J) 120 QR(K,KP1)=1.0+QR(K,KP1) GO TO 10 C C QR = ( F(X) - X + A, LAMBDA*DF(X) + (1 - LAMBDA)*I ) . C 140 DO 150 J=1,N 150 QR(J,1)=TZ(J)-Y(J+1)+A(J) DO 170 K=1,N CALL FJAC(Y(2),TZ,K) KP1=K+1 DO 160 J=1,N 160 QR(J,KP1)=Y(1)*TZ(J) 170 QR(K,KP1)=1.0-Y(1)+QR(K,KP1) C C * * * * * * * * * * * * * * * * * C REDUCE THE JACOBIAN MATRIX TO UPPER TRIANGULAR FORM. 10 CALL DCPOSE(NDIM,N,QR,ALPHA,PIVOT,IERR,TZ,YP) IF (IERR .EQ. 0) GO TO 20 IFLAG=4 RETURN C COMPUTE KERNAL OF JACOBIAN, WHICH SPECIFIES YP=DY/DS. 20 TZ(NP1)=1.0 DO 40 LW=1,N I=NP1-LW IK=I+1 SUM=0.0 DO 30 J=IK,NP1 30 SUM=SUM+QR(I,J)*TZ(J) 40 TZ(I)=-SUM/ALPHA(I) YPNORM=SQRT(INNER(1,NP1,TZ,TZ,0.0)) DO 60 K=1,NP1 KPIV=PIVOT(K) 60 YP(KPIV)=TZ(K)/YPNORM IF (INNER(1,NP1,YP,YPOLD,0.0) .GE. 0.0) GO TO 80 DO 70 I=1,NP1 70 YP(I)=-YP(I) C C SAVE CURRENT DERIVATIVE (= TANGENT VECTOR) IN YPOLD . 80 DO 90 I=1,NP1 90 YPOLD(I)=YP(I) RETURN END REAL FUNCTION INNER(K,M,A,B,C) DIMENSION A(1),B(1) SUM=C DO 10 I=K,M 10 SUM=SUM+A(I)*B(I) INNER=SUM RETURN END SUBROUTINE DCPOSE(NDIM,N,QR,ALPHA,PIVOT,IERR,Y,SUM) C C SUBROUTINE DCPOSE IS A MODIFICATION OF THE ALGOL PROCEDURE C DECOMPOSE IN P. BUSINGER AND G. H. GOLUB, LINEAR LEAST C SQUARES SOLUTIONS BY HOUSEHOLDER TRANSFORMATIONS, C NUMER. MATH. 7 (1965) 269-276. C INTEGER NDIM,N,PIVOT(1) REAL QR(NDIM,1),ALPHA(N) INTEGER IERR,I,J,JBAR,K REAL BETA,SIGMA,ALPHAK,QRKK,Y(1),SUM(1) REAL INNER IERR=0 NP1=N+1 DO 20 J=1,NP1 SUM(J)=INNER(1,N,QR(1,J),QR(1,J),0.0) 20 PIVOT(J)=J DO 500 K=1,N SIGMA=SUM(K) JBAR=K KP1=K+1 DO 40 J=KP1,NP1 IF (SIGMA .GE. SUM(J)) GO TO 40 SIGMA=SUM(J) JBAR=J 40 CONTINUE IF (JBAR .EQ. K) GO TO 70 I=PIVOT(K) PIVOT(K)=PIVOT(JBAR) PIVOT(JBAR)=I SUM(JBAR)=SUM(K) SUM(K)=SIGMA DO 50 I=1,N SIGMA=QR(I,K) QR(I,K)=QR(I,JBAR) QR(I,JBAR)=SIGMA 50 CONTINUE C END OF COLUMN INTERCHANGE. 70 SIGMA=INNER(K,N,QR(1,K),QR(1,K),0.0) IF (SIGMA .NE. 0.0) GO TO 60 IERR=1 RETURN 60 IF (K .EQ. N) GO TO 500 QRKK=QR(K,K) ALPHAK=-SQRT(SIGMA) IF (QRKK .LT. 0.0) ALPHAK=-ALPHAK ALPHA(K)=ALPHAK BETA=1.0/(SIGMA-QRKK*ALPHAK) QR(K,K)=QRKK-ALPHAK DO 80 J=KP1,NP1 80 Y(J)=BETA*INNER(K,N,QR(1,K),QR(1,J),0.0) DO 100 J=KP1,NP1 DO 90 I=K,N QR(I,J)=QR(I,J)-QR(I,K)*Y(J) 90 CONTINUE SUM(J)=SUM(J)-QR(K,J)**2 100 CONTINUE 500 CONTINUE ALPHA(N)=QR(N,N) RETURN END SUBROUTINE STEP(X,Y,F,NEQN,H,EPS,WT,START, 1 HOLD,K,KOLD,CRASH,PHI,P,YP,PSI) C C SUBROUTINE STEP INTEGRATES A SYSTEM OF FIRST ORDER ORDINARY C DIFFERENTIAL EQUATIONS ONE STEP NORMALLY FROM X TO X+H, USING A C MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS PECE FORMULAS. LOCAL C EXTRAPOLATION IS USED TO IMPROVE ABSOLUTE STABILITY AND ACCURACY. C THE CODE ADJUSTS ITS ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR C PER UNIT STEP IN A GENERALIZED SENSE. SPECIAL DEVICES ARE INCLUDED C TO CONTROL ROUNDOFF ERROR AND TO DETECT WHEN THE USER IS REQUESTING C TOO MUCH ACCURACY. C C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS& THE INITIAL C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON. C C C THE PARAMETERS REPRESENT& C X -- INDEPENDENT VARIABLE C Y(*) -- SOLUTION VECTOR AT X C YP(*) -- DERIVATIVE OF SOLUTION VECTOR AT X AFTER SUCCESSFUL C STEP C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C H -- APPROPRIATE STEP SIZE FOR NEXT STEP. NORMALLY DETERMINED BY C CODE C EPS -- LOCAL ERROR TOLERANCE. MUST BE VARIABLE C WT(*) -- VECTOR OF WEIGHTS FOR ERROR CRITERION C START -- LOGICAL VARIABLE SET .TRUE. FOR FIRST STEP, .FALSE. C OTHERWISE C HOLD -- STEP SIZE USED FOR LAST SUCCESSFUL STEP C K -- APPROPRIATE ORDER FOR NEXT STEP (DETERMINED BY CODE) C KOLD -- ORDER USED FOR LAST SUCCESSFUL STEP C CRASH -- LOGICAL VARIABLE SET .TRUE. WHEN NO STEP CAN BE TAKEN, C .FALSE. OTHERWISE C THE ARRAYS PHI, PSI ARE REQUIRED FOR THE INTERPOLATION SUBROUTINE C INTRP. THE ARRAY P IS INTERNAL TO THE CODE. C C INPUT TO STEP C C FIRST CALL -- C C THE USER MUST PROVIDE STORAGE IN HIS DRIVER PROGRAM FOR ALL ARRAYS C IN THE CALL LIST, NAMELY C C DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12) C C THE USER MUST ALSO DECLARE START AND CRASH LOGICAL VARIABLES C AND F AN EXTERNAL SUBROUTINE, SUPPLY THE SUBROUTINE F(X,Y,YP) C TO EVALUATE C DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN)) C AND INITIALIZE ONLY THE FOLLOWING PARAMETERS& C X -- INITIAL VALUE OF THE INDEPENDENT VARIABLE C Y(*) -- VECTOR OF INITIAL VALUES OF DEPENDENT VARIABLES C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED C H -- NOMINAL STEP SIZE INDICATING DIRECTION OF INTEGRATION C AND MAXIMUM SIZE OF STEP. MUST BE VARIABLE C EPS -- LOCAL ERROR TOLERANCE PER STEP. MUST BE VARIABLE C WT(*) -- VECTOR OF NON-ZERO WEIGHTS FOR ERROR CRITERION C START -- .TRUE. C C STEP REQUIRES THE L2 NORM OF THE VECTOR WITH COMPONENTS C LOCAL ERROR(L)/WT(L) BE LESS THAN EPS FOR A SUCCESSFUL STEP. THE C ARRAY WT ALLOWS THE USER TO SPECIFY AN ERROR TEST APPROPRIATE C FOR HIS PROBLEM. FOR EXAMPLE, C WT(L) = 1.0 SPECIFIES ABSOLUTE ERROR, C =ABS(Y(L)) ERROR RELATIVE TO THE MOST RECENT VALUE OF THE C L-TH COMPONENT OF THE DERIVATIVE, C = AMAX1(WT(L),ABS(Y(L))) ERROR RELATIVE TO THE LARGEST C MAGNITUDE OF L-TH COMPONENT OBTAINED SO FAR, C = ABS(Y(L))*RELERR/EPS + ABSERR/EPS SPECIFIES A MIXED C RELATIVE-ABSOLUTE TEST WHERE RELERR IS RELATIVE C ERROR, ABSERR IS ABSOLUTE ERROR AND EPS = C AMAX1(RELERR,ABSERR) C C SUBSEQUENT CALLS -- C C SUBROUTINE STEP IS DESIGNED SO THAT ALL INFORMATION NEEDED TO C CONTINUE THE INTEGRATION, INCLUDING THE STEP SIZE H AND THE ORDER C K, IS RETURNED WITH EACH STEP. WITH THE EXCEPTION OF THE STEP C SIZE, THE ERROR TOLERANCE, AND THE WEIGHTS, NONE OF THE PARAMETERS C SHOULD BE ALTERED. THE ARRAY WT MUST BE UPDATED AFTER EACH STEP C TO MAINTAIN RELATIVE ERROR TESTS LIKE THOSE ABOVE. NORMALLY THE C INTEGRATION IS CONTINUED JUST BEYOND THE DESIRED ENDPOINT AND THE C SOLUTION INTERPOLATED THERE WITH SUBROUTINE INTRP. IF IT IS C IMPOSSIBLE TO INTEGRATE BEYOND THE ENDPOINT, THE STEP SIZE MAY BE C REDUCED TO HIT THE ENDPOINT SINCE THE CODE WILL NOT TAKE A STEP C LARGER THAN THE H INPUT. CHANGING THE DIRECTION OF INTEGRATION, C I.E., THE SIGN OF H, REQUIRES THE USER SET START = .TRUE. BEFORE C CALLING STEP AGAIN. THIS IS THE ONLY SITUATION IN WHICH START C SHOULD BE ALTERED. C C OUTPUT FROM STEP C C SUCCESSFUL STEP -- C C THE SUBROUTINE RETURNS AFTER EACH SUCCESSFUL STEP WITH START AND C CRASH SET .FALSE. . X REPRESENTS THE INDEPENDENT VARIABLE C ADVANCED ONE STEP OF LENGTH HOLD FROM ITS VALUE ON INPUT AND Y C THE SOLUTION VECTOR AT THE NEW VALUE OF X . ALL OTHER PARAMETERS C REPRESENT INFORMATION CORRESPONDING TO THE NEW X NEEDED TO C CONTINUE THE INTEGRATION. C C UNSUCCESSFUL STEP -- C C WHEN THE ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION, C THE SUBROUTINE RETURNS WITHOUT TAKING A STEP AND CRASH = .TRUE. C AN APPROPRIATE STEP SIZE AND ERROR TOLERANCE FOR CONTINUING ARE C ESTIMATED AND ALL OTHER INFORMATION IS RESTORED AS UPON INPUT C BEFORE RETURNING. TO CONTINUE WITH THE LARGER TOLERANCE, THE USER C JUST CALLS THE CODE AGAIN. A RESTART IS NEITHER REQUIRED NOR C DESIRABLE. C LOGICAL START,CRASH,PHASE1,NORND DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12) DIMENSION ALPHA(12),BETA(12),SIG(13),W(12),V(12),G(13), 1 GSTR(13),TWO(13) EXTERNAL F C*********************************************************************** C* THE ONLY MACHINE DEPENDENT CONSTANTS ARE BASED ON THE MACHINE UNIT * C* ROUNDOFF ERROR U WHICH IS THE SMALLEST POSITIVE NUMBER SUCH THAT * C* 1.0+U .GT. 1.0 . THE USER MUST CALCULATE U AND INSERT * C* TWOU=2.0*U AND FOURU=4.0*U IN THE DATA STATEMENT BEFORE CALLING * C* THE CODE. THE ROUTINE MACHIN CALCULATES U. C* FOR THE CDC 6500, 6600, 7600 -- * C* DATA TWOU,FOURU /16424000000000000000B,16434000000000000000B/ C* FOR THE IBM 360/370 (SINGLE PRECISION) -- * C* DATA TWOU,FOURU /Z3C200000,Z3C400000/ C* FOR THE IBM 360/370 (DOUBLE PRECISION) -- * C* DATA TWOU,FOURU /Z3420000000000000,Z3440000000000000/ C* FOR THE UNIVAC 1108 AND HONEYWELL 6000 -- * C* DATA TWOU,FOURU /2.0**(-25),2.0**(-24)/ C* FOR THE PDP-11 -- * C* DATA TWOU,FOURU /2.0**(-22),2.0**(-21)/ C*********************************************************************** DATA TWO/2.0,4.0,8.0,16.0,32.0,64.0,128.0,256.0,512.0,1024.0, 1 2048.0,4096.0,8192.0/ DATA GSTR/0.500,0.0833,0.0417,0.0264,0.0188,0.0143,0.0114,0.00936, 1 0.00789,0.00679,0.00592,0.00524,0.00468/ DATA G(1),G(2)/1.0,0.5/,SIG(1)/1.0/ C C C *** BEGIN BLOCK 0 C CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE C PRECISION. IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A C STARTING STEP SIZE. C *** C C IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE C CRASH = .TRUE. IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 5 H = SIGN(FOURU*ABS(X),H) RETURN 5 P5EPS = 0.5*EPS C C IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE C ROUND = 0.0 DO 10 L = 1,NEQN 10 ROUND = ROUND + (Y(L)/WT(L))**2 ROUND = TWOU*SQRT(ROUND) IF(P5EPS .GE. ROUND) GO TO 15 EPS = 2.0*ROUND*(1.0 + FOURU) RETURN 15 CRASH = .FALSE. IF(.NOT.START) GO TO 99 C C INITIALIZE. COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP C CALL F(X,Y,YP) SUM = 0.0 DO 20 L = 1,NEQN PHI(L,1) = YP(L) PHI(L,2) = 0.0 20 SUM = SUM + (YP(L)/WT(L))**2 SUM = SQRT(SUM) ABSH = ABS(H) IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM) H = SIGN(AMAX1(ABSH,FOURU*ABS(X)),H) HOLD = 0.0 K = 1 KOLD = 0 START = .FALSE. PHASE1 = .TRUE. NORND = .TRUE. IF(P5EPS .GT. 100.0*ROUND) GO TO 99 NORND = .FALSE. DO 25 L = 1,NEQN 25 PHI(L,15) = 0.0 99 IFAIL = 0 C *** END BLOCK 0 *** C C *** BEGIN BLOCK 1 *** C COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP. AVOID COMPUTING C THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED. C *** C 100 KP1 = K+1 KP2 = K+2 KM1 = K-1 KM2 = K-2 C C NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT C ONE. WHEN K.LT.NS, NO COEFFICIENTS CHANGE IF(H .NE. HOLD) NS = 0 NS = MIN0(NS+1,KOLD+1) NSP1 = NS+1 IF (K .LT. NS) GO TO 199 C C COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH C ARE CHANGED C BETA(NS) = 1.0 REALNS = NS ALPHA(NS) = 1.0/REALNS TEMP1 = H*REALNS SIG(NSP1) = 1.0 IF(K .LT. NSP1) GO TO 110 DO 105 I = NSP1,K IM1 = I-1 TEMP2 = PSI(IM1) PSI(IM1) = TEMP1 BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2 TEMP1 = TEMP2 + H ALPHA(I) = H/TEMP1 REALI = I 105 SIG(I+1) = REALI*ALPHA(I)*SIG(I) 110 PSI(K) = TEMP1 C C COMPUTE COEFFICIENTS G(*) C C INITIALIZE V(*) AND SET W(*). G(2) IS SET IN DATA STATEMENT C IF(NS .GT. 1) GO TO 120 DO 115 IQ = 1,K TEMP3 = IQ*(IQ+1) V(IQ) = 1.0/TEMP3 115 W(IQ) = V(IQ) GO TO 140 C C IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*) C 120 IF(K .LE. KOLD) GO TO 130 TEMP4 = K*KP1 V(K) = 1.0/TEMP4 NSM2 = NS-2 IF(NSM2 .LT. 1) GO TO 130 DO 125 J = 1,NSM2 I = K-J 125 V(I) = V(I) - ALPHA(J+1)*V(I+1) C C UPDATE V(*) AND SET W(*) C 130 LIMIT1 = KP1 - NS TEMP5 = ALPHA(NS) DO 135 IQ = 1,LIMIT1 V(IQ) = V(IQ) - TEMP5*V(IQ+1) 135 W(IQ) = V(IQ) G(NSP1) = W(1) C C COMPUTE THE G(*) IN THE WORK VECTOR W(*) C 140 NSP2 = NS + 2 IF(KP1 .LT. NSP2) GO TO 199 DO 150 I = NSP2,KP1 LIMIT2 = KP2 - I TEMP6 = ALPHA(I-1) DO 145 IQ = 1,LIMIT2 145 W(IQ) = W(IQ) - TEMP6*W(IQ+1) 150 G(I) = W(1) 199 CONTINUE C *** END BLOCK 1 *** C C *** BEGIN BLOCK 2 *** C PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED C SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K, C K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED. C *** C C CHANGE PHI TO PHI STAR C IF(K .LT. NSP1) GO TO 215 DO 210 I = NSP1,K TEMP1 = BETA(I) DO 205 L = 1,NEQN 205 PHI(L,I) = TEMP1*PHI(L,I) 210 CONTINUE C C PREDICT SOLUTION AND DIFFERENCES C 215 DO 220 L = 1,NEQN PHI(L,KP2) = PHI(L,KP1) PHI(L,KP1) = 0.0 220 P(L) = 0.0 DO 230 J = 1,K I = KP1 - J IP1 = I+1 TEMP2 = G(I) DO 225 L = 1,NEQN P(L) = P(L) + TEMP2*PHI(L,I) 225 PHI(L,I) = PHI(L,I) + PHI(L,IP1) 230 CONTINUE IF(NORND) GO TO 240 DO 235 L = 1,NEQN TAU = H*P(L) - PHI(L,15) P(L) = Y(L) + TAU 235 PHI(L,16) = (P(L) - Y(L)) -TAU GO TO 250 240 DO 245 L = 1,NEQN 245 P(L) = Y(L) + H*P(L) 250 XOLD = X X = X + H ABSH = ABS(H) CALL F(X,P,YP) C C ESTIMATE ERRORS AT ORDERS K,K-1,K-2 ERKM2 = 0.0 ERKM1 = 0.0 ERK = 0.0 DO 265 L = 1,NEQN TEMP3 = 1.0/WT(L) TEMP4 = YP(L) - PHI(L,1) IF(KM2) 265,260,255 255 ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2 260 ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2 265 ERK = ERK + (TEMP4*TEMP3)**2 IF(KM2) 280,275,270 270 ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*SQRT(ERKM2) 275 ERKM1 = ABSH*SIG(K)*GSTR(KM1)*SQRT(ERKM1) 280 TEMP5 = ABSH*SQRT(ERK) ERR = TEMP5*(G(K)-G(KP1)) ERK = TEMP5*SIG(KP1)*GSTR(K) KNEW = K C C TEST IF ORDER SHOULD BE LOWERED C IF(KM2) 299,290,285 285 IF(AMAX1(ERKM1,ERKM2) .LE. ERK) KNEW = KM1 GO TO 299 290 IF(ERKM1 .LE. 0.5*ERK) KNEW = KM1 C C TEST IF STEP SUCCESSFUL C 299 IF(ERR .LE. EPS) GO TO 400 C *** END BLOCK 2 *** C C *** BEGIN BLOCK 3 *** C THE STEP IS UNSUCCESSFUL. RESTORE X, PHI(*,*), PSI(*) C IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE. IF STEP FAILS MORE C THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE. DOUBLE ERROR C TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE C PRECISION C *** C C RESTORE X, PHI(*,*) AND PSI(*) C PHASE1 = .FALSE. X = XOLD DO 310 I = 1,K TEMP1 = 1.0/BETA(I) IP1 = I+1 DO 305 L = 1,NEQN 305 PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1)) 310 CONTINUE IF(K .LT. 2) GO TO 320 DO 315 I = 2,K 315 PSI(I-1) = PSI(I) - H C C ON THIRD FAILURE, SET ORDER TO ONE. THEREAFTER, USE OPTIMAL STEP C SIZE C 320 IFAIL = IFAIL + 1 TEMP2 = 0.5 IF(IFAIL - 3) 335,330,325 325 IF(P5EPS .LT. 0.25*ERK) TEMP2 = SQRT(P5EPS/ERK) 330 KNEW = 1 335 H = TEMP2*H K = KNEW IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 340 CRASH = .TRUE. H = SIGN(FOURU*ABS(X),H) EPS = EPS + EPS RETURN 340 GO TO 100 C *** END BLOCK 3 *** C C *** BEGIN BLOCK 4 *** C THE STEP IS SUCCESSFUL. CORRECT THE PREDICTED SOLUTION, EVALUATE C THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE C DIFFERENCES. DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP. C *** 400 KOLD = K HOLD = H C C CORRECT AND EVALUATE C TEMP1 = H*G(KP1) IF(NORND) GO TO 410 DO 405 L = 1,NEQN RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16) Y(L) = P(L) + RHO 405 PHI(L,15) = (Y(L) - P(L)) - RHO GO TO 420 410 DO 415 L = 1,NEQN 415 Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1)) 420 CALL F(X,Y,YP) C C UPDATE DIFFERENCES FOR NEXT STEP C DO 425 L = 1,NEQN PHI(L,KP1) = YP(L) - PHI(L,1) 425 PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2) DO 435 I = 1,K DO 430 L = 1,NEQN 430 PHI(L,I) = PHI(L,I) + PHI(L,KP1) 435 CONTINUE C C ESTIMATE ERROR AT ORDER K+1 UNLESS& C IN FIRST PHASE WHEN ALWAYS RAISE ORDER, C ALREADY DECIDED TO LOWER ORDER, C STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE C ERKP1 = 0.0 IF(KNEW .EQ. KM1 .OR. K .EQ. 12) PHASE1 = .FALSE. IF(PHASE1) GO TO 450 IF(KNEW .EQ. KM1) GO TO 455 IF(KP1 .GT. NS) GO TO 460 DO 440 L = 1,NEQN 440 ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2 ERKP1 = ABSH*GSTR(KP1)*SQRT(ERKP1) C C USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER C FOR NEXT STEP C IF(K .GT. 1) GO TO 445 IF(ERKP1 .GE. 0.5*ERK) GO TO 460 GO TO 450 445 IF(ERKM1 .LE. AMIN1(ERK,ERKP1)) GO TO 455 IF(ERKP1 .GE. ERK .OR. K .EQ. 12) GO TO 460 C C HERE ERKP1 .LT. ERK .LT. AMAX1(ERKM1,ERKM2) ELSE ORDER WOULD HAVE C BEEN LOWERED IN BLOCK 2. THUS ORDER IS TO BE RAISED C C RAISE ORDER C 450 K = KP1 ERK = ERKP1 GO TO 460 C C LOWER ORDER C 455 K = KM1 ERK = ERKM1 C C WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP C 460 HNEW = H + H IF(PHASE1) GO TO 465 IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465 HNEW = H IF(P5EPS .GE. ERK) GO TO 465 TEMP2 = K+1 R = (P5EPS/ERK)**(1.0/TEMP2) HNEW = ABSH*AMAX1(0.5,AMIN1(0.9,R)) HNEW = SIGN(AMAX1(HNEW,FOURU*ABS(X)),H) 465 H = HNEW RETURN C *** END BLOCK 4 *** END SUBROUTINE INTRP(X,Y,XOUT,YOUT,YPOUT,NEQN,KOLD,PHI,PSI) C C THE METHODS IN SUBROUTINE STEP APPROXIMATE THE SOLUTION NEAR X C BY A POLYNOMIAL. SUBROUTINE INTRP APPROXIMATES THE SOLUTION AT C XOUT BY EVALUATING THE POLYNOMIAL THERE. INFORMATION DEFINING THIS C POLYNOMIAL IS PASSED FROM STEP SO INTRP CANNOT BE USED ALONE. C C THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT, C COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS& THE INITIAL C VALUE PROBLEM BY L. F. SHAMPINE AND M. K. GORDON. C C INPUT TO INTRP -- C C THE USER PROVIDES STORAGE IN THE CALLING PROGRAM FOR THE ARRAYS IN C THE CALL LIST DIMENSION Y(NEQN),YOUT(NEQN),YPOUT(NEQN),PHI(NEQN,16),PSI(12) C AND DEFINES C XOUT -- POINT AT WHICH SOLUTION IS DESIRED. C THE REMAINING PARAMETERS ARE DEFINED IN STEP AND PASSED TO INTRP C FROM THAT SUBROUTINE C C OUTPUT FROM INTRP -- C C YOUT(*) -- SOLUTION AT XOUT C YPOUT(*) -- DERIVATIVE OF SOLUTION AT XOUT C THE REMAINING PARAMETERS ARE RETURNED UNALTERED FROM THEIR INPUT C VALUES. INTEGRATION WITH STEP MAY BE CONTINUED. C DIMENSION G(13),W(13),RHO(13) DATA G(1)/1.0/,RHO(1)/1.0/ C HI = XOUT - X KI = KOLD + 1 KIP1 = KI + 1 C C INITIALIZE W(*) FOR COMPUTING G(*) C DO 5 I = 1,KI TEMP1 = I 5 W(I) = 1.0/TEMP1 TERM = 0.0 C C COMPUTE G(*) C DO 15 J = 2,KI JM1 = J - 1 PSIJM1 = PSI(JM1) GAMMA = (HI + TERM)/PSIJM1 ETA = HI/PSIJM1 LIMIT1 = KIP1 - J DO 10 I = 1,LIMIT1 10 W(I) = GAMMA*W(I) - ETA*W(I+1) G(J) = W(1) RHO(J) = GAMMA*RHO(JM1) 15 TERM = PSIJM1 C C INTERPOLATE C DO 20 L = 1,NEQN YPOUT(L) = 0.0 20 YOUT(L) = 0.0 DO 30 J = 1,KI I = KIP1 - J TEMP2 = G(I) TEMP3 = RHO(I) DO 25 L = 1,NEQN YOUT(L) = YOUT(L) + TEMP2*PHI(L,I) 25 YPOUT(L) = YPOUT(L) + TEMP3*PHI(L,I) 30 CONTINUE DO 35 L = 1,NEQN 35 YOUT(L) = Y(L) + HI*YOUT(L) RETURN END SUBROUTINE ROOT(T,FT,B,C,RELERR,ABSERR,IFLAG) C C ROOT COMPUTES A ROOT OF THE NONLINEAR EQUATION F(X)=0 C WHERE F(X) IS A CONTINOUS REAL FUNCTION OF A SINGLE REAL C VARIABLE X. THE METHOD USED IS A COMBINATION OF BISECTION C AND THE SECANT RULE. C C NORMAL INPUT CONSISTS OF A CONTINUOS FUNCTION F AND AN C INTERVAL (B,C) SUCH THAT F(B)*F(C).LE.0.0. EACH ITERATION C FINDS NEW VALUES OF B AND C SUCH THAT THE INTERVAL(B,C) IS C SHRUNK AND F(B)*F(C).LE.0.0. THE STOPPING CRITERION IS C C ABS(B-C).LE.2.0*(RELERR*ABS(B)+ABSERR) C C WHERE RELERR=RELATIVE ERROR AND ABSERR=ABSOLUTE ERROR ARE C INPUT QUANTITIES. SET THE FLAG, IFLAG, POSITIVE TO INITIALIZE C THE COMPUTATION. AS B,C AND IFLAG ARE USED FOR BOTH INPUT AND C OUTPUT, THEY MUST BE VARIABLES IN THE CALLING PROGRAM. C C IF 0 IS A POSSIBLE ROOT, ONE SHOULD NOT CHOOSE ABSERR=0.0. C C THE OUTPUT VALUE OF B IS THE BETTER APPROXIMATION TO A ROOT C AS B AND C ARE ALWAYS REDEFINED SO THAT ABS(F(B)).LE.ABS(F(C)). C C TO SOLVE THE EQUATION, ROOT MUST EVALUATE F(X) REPEATEDLY. THIS C IS DONE IN THE CALLING PROGRAM. WHEN AN EVALUATION OF F IS C NEEDED AT T, ROOT RETURNS WITH IFLAG NEGATIVE. EVALUATE FT=F(T) C AND CALL ROOT AGAIN. DO NOT ALTER IFLAG. C C WHEN THE COMPUTATION IS COMPLETE, ROOT RETURNS TO THE CALLING C PROGRAM WITH IFLAG POSITIVE= C C IFLAG=1 IF F(B)*F(C).LT.0 AND THE STOPPING CRITERION IS MET. C C =2 IF A VALUE B IS FOUND SUCH THAT THE COMPUTED VALUE C F(B) IS EXACTLY ZERO. THE INTERVAL (B,C) MAY NOT C SATISFY THE STOPPING CRITERION. C C =3 IF ABS(F(B)) EXCEEDS THE INPUT VALUES ABS(F(B)), C ABS(F(C)). IN THIS CASE IT IS LIKELY THAT B IS CLOSE C TO A POLE OF F. C C =4 IF NO ODD ORDER ROOT WAS FOUND IN THE INTERVAL. A C LOCAL MINIMUM MAY HAVE BEEN OBTAINED. C C =5 IF TOO MANY FUNCTION EVALUATIONS WERE MADE. C (AS PROGRAMMED, 500 ARE ALLOWED.) C C THIS CODE IS A MODIFICATION OF THE CODE ZEROIN WHICH IS COMPLETELY C EXPLAINED AND DOCUMENTED IN THE TEXT, NUMERICAL COMPUTING= AN C INTRODUCTION BY L. F. SHAMPINE AND R. C. ALLEN. C C*********************************************************************** C* THE ONLY MACHINE DEPENDENT CONSTANT IS BASED ON THE MACHINE UNIT * C* ROUNDOFF ERROR U WHICH IS THE SMALLEST POSITIVE NUMBER SUCH THAT * C* 1.0+U .GT. 1.0. U MUST BE CALCULATED AND INSERTED IN THE * C* FOLLOWING DATA STATEMENT BEFORE USING ROOT. THE ROUTINE MACHIN * C* CALCULATES U. C* FOR THE CDC 6500, 6600, 7600 -- * C* DATA U /16414000000000000000B/ C* FOR THE IBM 360/370 (SINGLE PRECISION) -- * C* DATA U /Z3C100000/ C* FOR THE IBM 360/370 (DOUBLE PRECISION) -- * C* DATA U /Z3410000000000000/ C* FOR THE UNIVAC 1108 AND HONEYWELL 6000 -- * C* DATA U /2.0**(-26)/ C* FOR THE PDP-11 -- * C* DATA U /2.0**(-23)/ C*********************************************************************** C IF(IFLAG.GE.0) GO TO 100 IFLAG=IABS(IFLAG) GO TO (200,300,400), IFLAG 100 RE=AMAX1(RELERR,U) AE=AMAX1(ABSERR,0.0) IC=0 ACBS=ABS(B-C) A=C T=A IFLAG=-1 RETURN 200 FA=FT T=B IFLAG=-2 RETURN 300 FB=FT FC=FA KOUNT=2 FX=AMAX1(ABS(FB),ABS(FC)) 1 IF(ABS(FC).GE.ABS(FB))GO TO 2 C C INTERCHANGE B AND C SO YHAT ABS(F(B)).LE.ABS(F(C)). C A=B FA=FB B=C FB=FC C=A FC=FA 2 CMB=0.5*(C-B) ACMB=ABS(CMB) TOL=RE*ABS(B)+AE C C TEST STOPPING CRITERION AND FUNCTION COUNT. C IF(ACMB.LE.TOL)GO TO 8 IF(KOUNT.GE.500)GO TO 12 C C CALCULATE NEW ITERATE EXPLICITLY AS B+P/Q C WHERE WE ARRANGE P.GE.0. THE IMPLICIT C FORM IS USED TO PREVENT OVERFLOW. C P=(B-A)*FB Q=FA-FB IF(P.GE.0.0)GO TO 3 P=-P Q=-Q C C UPDATE A, CHECK IF REDUCTION IN THE SIZE OF BRACKETING C INTERVAL IS SATISFACTORY. IF NOT BISECT UNTIL IT IS. C 3 A=B FA=FB IC=IC+1 IF(IC.LT.4)GO TO 4 IF(8.0*ACMB.GE.ACBS)GO TO 6 IC=0 ACBS=ACMB C C TEST FOR TOO SMALL A CHANGE. C 4 IF(P.GT.ABS(Q)*TOL)GO TO 5 C C INCREMENT BY TOLERANCE C B=B+SIGN(TOL,CMB) GO TO 7 C C ROOT OUGHT TO BE BETWEEN B AND (C+B)/2 C 5 IF(P.GE.CMB*Q)GO TO 6 C C USE SECANT RULE. C B=B+P/Q GO TO 7 C C USE BISECTION. C 6 B=0.5*(C+B) C C HAVE COMPLETED COMPUTATION FOR NEW ITERATE B. C 7 T=B IFLAG=-3 RETURN 400 FB=FT IF(FB.EQ.0.0)GO TO 9 KOUNT=KOUNT+1 IF(SIGN(1.0,FB).NE.SIGN(1.0,FC))GO TO 1 C=A FC=FA GO TO 1 C C FINISHED. SET IFLAG. C 8 IF(SIGN(1.0,FB).EQ.SIGN(1.0,FC))GO TO 11 IF(ABS(FB).GT.FX)GO TO 10 IFLAG=1 RETURN 9 IFLAG=2 RETURN 10 IFLAG=3 RETURN 11 IFLAG=4 RETURN 12 IFLAG=5 RETURN END