SUBROUTINE RKF45(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK,IWORK, * CUTOFF,KT) Implicit DoublePrecision (A-H,O-Z) CCCCC CCCCCRKF45 WITH CUTOFF FEATURE. SEE RUNCUT.DOC FOR DOCUMENTATION. CCCCCIN THIS VERSION, 3 NEW PARAMETERS HAVE BEEN ADDED TO THE CCCCCINTEGER ARRAY IWORK. THEY ARE... CCCCCNSUCC = IWORK(6) = NUMBER OF SUCCESSFUL STEPS CCCCCNFAIL = IWORK(7) = NUMBER OF FAILED STEP ATTEMPTS CCCCCNSWF = IWORK(8) = NUMBER OF STEPS ON WHICH AT LEAST CCCCC ONE FAILURE OCCURED BEFORE SUCCESS CCCCCTHESE ARE ALL INITIALIZED TO ZERO ON THE FIRST CALL TO RKF45 CCCCCBUT NOT ON SUBSEQUENT CALLS, SO THEY CONTINUE TO ACCUMULATE. CCCCCTHE TOTAL NUMBER OF FUNCTION EVALUATIONS = NFE = IWORK(1) IS CCCCC RELATED (APPROXIMATELY) BY: NFE = NSUCC*6 + NFAIL*5 SINCE CCCCC A SUCCESSFUL STEP COSTS 6 EVALUATIONS AND A FAILED STEP ATTEMPT CCCCC COSTS 5 EVALUATIONS. CCCCC CCCCCON RETURN, THE VALUE OF THE STEP SIZE THAT WILL BE ATTEMPTED CCCCC ON THE NEXT STEP IS CONTAINED IN WORK(NEQN+1). THIS IS CCCCCA STANDARD PART OF THE ORIGINAL RKF45 CODE. IN THIS VERSION, CCCCC ONE MAY ALSO RETRIEVE THE VALUE OF THE MOST RECENTLY USED CCCCC STEP SIZE, WHICH HAS BEEN PUT INTO WORK(NEQN+2). CCCCC INTEGER NEQN,IFLAG,IWORK(8) DoublePrecision Y(NEQN),T,TOUT,RELERR,ABSERR,WORK(1) C IF COMPILER CHECKS SUBSCRIPTS, CHANGE WORK(1) TO WORK(3+6*NEQN) C********** C NEW PARAMETERS DoublePrecision CUTOFF INTEGER KT C********** EXTERNAL F INTEGER K1,K2,K3,K4,K5,K6,K1M C C COMPUTE INDICES FOR THE SPLITTING OF THE WORK ARRAY C K1M=NEQN+1 K1=K1M+1 K2=K1+NEQN K3=K2+NEQN K4=K3+NEQN K5=K4+NEQN K6=K5+NEQN CCCCC CCCCC NOTICE 3 EXTRA PARAMETERS IWORK 6,7,8. CALL RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,WORK(1),WORK(K1M), 1 WORK(K1),WORK(K2),WORK(K3),WORK(K4),WORK(K5),WORK(K6), 2 WORK(K6+1),IWORK(1),IWORK(2),IWORK(3),IWORK(4),IWORK(5), * IWORK(6),IWORK(7),IWORK(8),CUTOFF,KT) RETURN END C C END OF INTERFACE ROUTINE. START INTEGRATION ROUTINE. C CCCCC CCCCC NOTICE EXTRA PARAMETERS... SUBROUTINE RKFS(F,NEQN,Y,T,TOUT,RELERR,ABSERR,IFLAG,YP,H,F1,F2,F3, 1 F4,F5,SAVRE,SAVAE,NFE,KOP,INIT,JFLAG,KFLAG, * NSUCC,NFAIL,NSWF,CUTOFF,KT) Implicit DoublePrecision (A-H,O-Z) C YP - DERIVATIVE OF SOLUTION VECTOR AT T C H - AN APPROPRIATE STEPSIZE TO BE USED FOR THE NEXT STEP C NFE - COUNTER ON THE NUMBER OF DERIVATIVE FUNCTION EVALUATIONS C LOGICAL HFAILD, OUTPUT INTEGER NEQN,IFLAG,NFE,KOP,INIT,JFLAG,KFLAG DoublePrecision Y(NEQN),T,TOUT,RELERR,ABSERR,H,YP(NEQN), 1 F1(NEQN),F2(NEQN),F3(NEQN),F4(NEQN),F5(NEQN),SAVRE,SAVAE C********** DoublePrecision CUTOFF INTEGER KT LOGICAL SKIP,CROSS CCCCC CCCCC NEW PARAMETERS... INTEGER NSUCC,NFAIL,NSWF C********** EXTERNAL F DoublePrecision A,AE,DT,EE,EEOET,ESTTOL,ET,HMIN,REMIN,RER,S, 1 SCALE,TOL,TOLN,U26,EPSP1,EPS,YPK INTEGER K,MAXNFE,MFLAG DoublePrecision dMAX1,dMIN1 DATA REMIN/1.E-12/ DATA MAXNFE/10000/ C********** C 'SKIP' DETERMINES WHETHER CUTOFF FEATURE IS ENABLED OR DISABLED. SKIP=((KT.LT.1).OR.(KT.GT.NEQN)) C********** C C CHECK INPUT PARAMETERS C IF(NEQN.LT.1) GO TO 10 IF((RELERR.LT.0.0) .OR. (ABSERR.LT.0.0)) GO TO 10 MFLAG=IABS(IFLAG) IF ((MFLAG.EQ. 0) .OR. (MFLAG.GT.8)) GO TO 10 IF(MFLAG.NE. 1) GO TO 20 C C FIRST CALL, COMPUTE MACHINE EPSILON C EPS=1.0 5 EPS=EPS/2.0 EPSP1=EPS+1.0 IF (EPSP1.GT.1.0) GO TO 5 U26=26.0*EPS GO TO 50 10 IFLAG=8 RETURN C C CHECK CONTINUATION POSSIBILITIES C 20 IF((T.EQ.TOUT) .AND. (KFLAG.NE.3)) GO TO 10 IF(MFLAG.NE. 2) GO TO 25 C IFLAG = +2 OR -2 IF ((KFLAG .EQ. 3) .OR.(INIT.EQ.0)) GO TO 45 IF (KFLAG.EQ.4) GO TO 40 IF ((KFLAG.EQ.5) .AND. (ABSERR.EQ.0.0)) GO TO 30 IF ((KFLAG.EQ.6) .AND. (RELERR.LE.SAVRE) .AND. 1 (ABSERR.LE.SAVAE)) GO TO 30 GO TO 50 C C IFLAG=3,4,5,6,7 OR 8 25 IF(IFLAG.EQ.3) GO TO 45 IF(IFLAG.EQ.4) GO TO 40 IF ((IFLAG.EQ.5) .AND. (ABSERR.GT.0.0)) GO TO 45 30 STOP 40 NFE=0 IF (MFLAG.EQ.2) GO TO 50 45 IFLAG=JFLAG IF(KFLAG.EQ.3) MFLAG=IABS(IFLAG) 50 JFLAG=IFLAG KFLAG=0 SAVRE=RELERR SAVAE=ABSERR RER=2.0*EPS+REMIN IF (RELERR.GE.RER) GO TO 55 RELERR=RER IFLAG=3 KFLAG=3 RETURN 55 DT=TOUT-T IF(MFLAG.EQ.1) GO TO 60 IF(INIT.EQ.0) GO TO 65 GO TO 80 C C INITIALIZATION C SET INITIALIZATION COMPLETION INDICATOR,INIT C SET INDICATOR FOR TOO MANY OUTPUT POINTS,KOP C EVALUATE INITIAL DERIVATIVES C SET COUNTER FOR FUNCTION EVALUATIONS, NFE C ESTIMATE STARTING STEPSIZE 60 INIT=0 KOP=0 CCCCC CCCCC INITIALIZE NEW PARAMETERS. CCCCC NOTICE THAT THEY ARE CUMULATIVE (LIKE NFE) CCCCC I.E. NOT RESET WITH EACH NEW CALL TO RKF45 NSUCC=0 NFAIL=0 NSWF =0 CCCCC A=T CALL F(A,Y,YP) NFE=1 IF(T.NE.TOUT)GO TO 65 IFLAG=2 RETURN 65 INIT=1 H=ABS(DT) TOLN=0. DO 70 K=1,NEQN TOL=RELERR*ABS(Y(K))+ABSERR IF(TOL.LE.0.)GO TO 70 TOLN=TOL YPK=ABS(YP(K)) IF(YPK*H**5.GT.TOL)H=(TOL/YPK)**0.2 70 CONTINUE IF(TOLN.LE.0.0)H=0.0 H=dMAX1(H,U26*dMAX1(ABS(T),ABS(DT))) JFLAG=ISIGN(2,IFLAG) 80 H=SIGN(H,DT) IF (ABS(H).GE.2.0*ABS(DT)) KOP=KOP+1 IF (KOP.NE.100) GO TO 85 KOP=0 IFLAG=7 RETURN 85 IF(ABS(DT).GT.U26*ABS(T)) GO TO 95 DO 90 K=1,NEQN 90 Y(K)=Y(K)+DT*YP(K) A=TOUT CALL F(A,Y,YP) NFE=NFE+1 GO TO 300 95 OUTPUT=.FALSE. SCALE=2.0/RELERR AE=SCALE*ABSERR 100 HFAILD=.FALSE. HMIN=U26*ABS(T) DT=TOUT-T IF(ABS(DT).GE.2.0*ABS(H)) GO TO 200 IF(ABS(DT).GT.ABS(H)) GO TO 150 OUTPUT=.TRUE. H=DT GO TO 200 150 H=0.5*DT 200 IF (NFE.LE.MAXNFE) GO TO 220 C C TOO MUCH WORK IFLAG=4 KFLAG=4 RETURN 220 CALL FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,F1) NFE=NFE+5 EEOET=0.0 DO 250 K=1,NEQN ET=ABS(Y(K))+ABS(F1(K))+AE IF(ET.GT.0.0) GO TO 240 C INAPPROPRIATE ERROR TOLERANCE IFLAG=5 RETURN 240 EE=ABS((-2090.0*YP(K)+(21970.0*F3(K)-15048.0*F4(K)))+ 1 (22528.0*F2(K)-27360.0*F5(K))) 250 EEOET=dMAX1(EEOET,EE/ET) ESTTOL=ABS(H)*EEOET*SCALE/7.524d5 IF(ESTTOL.LE.1.0) GO TO 260 C CUNSUCCESSFUL STEP C CCCCC CCCCC NFAIL COUNTS NUMBER OF FAILURES NFAIL=NFAIL+1 CCCCC HFAILD=.TRUE. OUTPUT=.FALSE. S=0.1 IF(ESTTOL.LT.5.9049d4)S=0.9/ESTTOL**0.2 H=S*H IF(ABS(H).GT.HMIN)GO TO 200 C C REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE IFLAG=6 KFLAG=6 RETURN 260 CONTINUE C CSUCCESSFUL STEP C CCCCC CCCCC NSUCC COUNTS NUMBER OF SUCCESSFUL STEPS. CCCCC NSWF INCREMENTS IF THIS SUCCESSFUL STEP WAS PRECEEDED BY CCCCC AT LEAST ONE STEP FAILURE. NSUCC=NSUCC+1 IF (HFAILD) NSWF=NSWF+1 CCCCC C********** C IF F1(KT) AND Y(KT) ARE ON OPPOSITE SIDES OF THE CUTOFF VALUE C THEN INTERPOLATE SOLUTION AND RETURN, ELSE CONTINUE INTEGRATING. IF (SKIP.OR..NOT.CROSS(Y(KT),YP(KT),F1(KT),CUTOFF)) GO TO 269 CALL INTERP(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,CUTOFF,KT,NFE) GO TO 272 269 CONTINUE C END OF INTERPOLATION SECTION. C********** T=T+H DO 270 K=1,NEQN 270 Y(K)=F1(K) A=T CALL F(A,Y,YP) NFE=NFE+1 C********** C SPECIAL CASE: IF Y(KT) HITS CUTOFF EXACTLY, SET T*=T, Y*=Y. IF (SKIP.OR.(Y(KT).NE.CUTOFF)) GO TO 274 DO 271 K=1,NEQN 271 F4(K)=Y(K) F3(NEQN)=T 272 IFLAG=ISIGN(9,IFLAG) 274 CONTINUE C********** CCCCC CCCCC SAVE STEP SIZE JUST USED IN F5(1) F5(1)=H CCCCC S=5.0 IF(ESTTOL.GT.1.889568E-4)S=0.9/ESTTOL**0.2 IF(HFAILD) S=dMIN1(S,1.d0) H=SIGN(dMAX1(S*ABS(H),HMIN),H) C********** CRETURN IF INTERPOLATION WAS DONE IF (ABS(IFLAG).EQ.9) RETURN C********** IF (OUTPUT) GO TO 300 IF(IFLAG.GT.0) GO TO 100 IFLAG=-2 RETURN 300 T=TOUT IFLAG=2 RETURN END C C ONE STEP INTEGRATOR. C SUBROUTINE FEHL(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,S) Implicit DoublePrecision (A-H,O-Z) INTEGER NEQN DoublePrecision Y(NEQN),T,H,YP(NEQN),F1(NEQN),F2(NEQN), 1 F3(NEQN),F4(NEQN),F5(NEQN),S(NEQN) DoublePrecision CH INTEGER K CH=H/4.0 DO 221 K=1,NEQN 221 F5(K)=Y(K)+CH*YP(K) CALL F(T+CH,F5,F1) CH=3.0*H/32.0 DO 222 K=1,NEQN 222 F5(K)=Y(K)+CH*(YP(K)+3.0*F1(K)) CALL F(T+3.0*H/8.0,F5,F2) CH=H/2197.0 DO 223 K=1,NEQN 223 F5(K)=Y(K)+CH*(1932.0*YP(K)+(7296*F2(K)-7200.0*F1(K))) CALL F(T+12.0*H/13.0,F5,F3) CH=H/4104.0 DO 224 K=1,NEQN 224 F5(K)=Y(K)+CH*((8341.0*YP(K)-845.0*F3(K))+ 1 (29440.0*F2(K)-3.2832d4*F1(K))) CALL F(T+H,F5,F4) CH=H/20520.0 DO 225 K=1,NEQN 225 F1(K)=Y(K)+CH*((-6080*YP(K)+(9295.0*F3(K)- 1 5643.0*F4(K)))+(4.1040d4*F1(K)-28352.0*F2(K))) CALL F(T+H/2.0,F1,F5) CH=H/7.618050d6 DO 230 K=1,NEQN 230 S(K)=Y(K)+CH*((9.02880d5*YP(K)+(3.855735d6*F3(K)- 1 1.371249d6*F4(K)))+(3.953664d6*F2(K)+2.77020d5*F5(K))) RETURN END C GENERAL BISECTION ROUTINE FOR SOLVING F(X) = 0. C ASSUMES ( AA .LT. BB ) .AND. ( F(AA) * F(BB) .LT. 0.0 ). C ACTUAL PRECISION = MAX ( TOL, MACHINE PRECISION ). DoublePrecision FUNCTION BI(AA,BB,F,TOL) Implicit DoublePrecision (A-H,O-Z) EXTERNAL F A=AA FA=F(A) BI=BB 10 B=BI 20 E=(B-A)/2.0 BI=A+E IF ((E.LE.TOL).OR.(BI.LE.A).OR.(BI.GE.B)) RETURN FM=F(BI) IF (FM.EQ.0.0) RETURN IF (FA*FM.LE.0.0) GO TO 10 A=BI FA=FM GO TO 20 END C C 'CROSS' RETURNS .TRUE. IF CUTOFF LEVEL IS CROSSED. LOGICAL FUNCTION CROSS(Y1,Y1P,Y2,CUT) Implicit DoublePrecision (A-H,O-Z) D1=Y1-CUT D2=Y2-CUT IF (D1) 10,20,30 10 IF (D2) 50,50,40 20 IF (Y1P) 10,50,30 30 IF (D2) 40,50,50 40 CROSS=.TRUE. RETURN 50 CROSS=.FALSE. RETURN END SUBROUTINE INTERP(F,NEQN,Y,T,H,YP,F1,F2,F3,F4,F5,CUTOFF,KT,NFE) Implicit DoublePrecision (A-H,O-Z) C CUBIC HERMITE INTERPOLATION INTEGER NEQN,KT,NFE,K DoublePrecision Y(NEQN),T,H,YP(NEQN),F1(NEQN),F2(NEQN), 1 F3(NEQN),F4(NEQN), F5(NEQN),CUTOFF,A,S COMMON /CFTS/ COEF(4),H8 EXTERNAL F,GETCO,BI,EVAL H8=H/8.0 C C F1 HAS SOLUTION VALUES AT T+H. GET SLOPES AT T+H. A=T+H CALL F(A,F1,F2) NFE=NFE+1 C C GET COEFFICIENTS OF THE INTERPOLANT OF THE 'KT' COMPONENT CALL GETCO(Y(KT),YP(KT),F1(KT),F2(KT)) COEF(1)=COEF(1)-CUTOFF C C SOLVE P(S) - CUTOFF = 0 BY BISECTION. NOTE -1 < S < 1. S=BI(-1.0,1.0,EVAL,0.0) C C INTERPOLATE OTHER COMPONENTS AND EVALUATE AT S DO 40 K=1,NEQN IF (K.EQ.KT) GO TO 40 CALL GETCO(Y(K),YP(K),F1(K),F2(K)) F4(K)=EVAL(S) 40 CONTINUE F4(KT)=CUTOFF C C COMPUTE T* FROM S BY REVERSING THE SCALING. F3(NEQN)=T+H*(S+1.0)/2.0 C C UPDATE Y AND YP TO NEW VALUES AT T+H T=A DO 50 K=1,NEQN Y(K) =F1(K) YP(K)=F2(K) 50 CONTINUE RETURN END C C THE CUBIC INTERPOLANT ON THE INTERVAL [T1,T2] IS C P(S) = C(1) + C(2) * S + C(3) * S**2 + C(4) * S**3 C WHERE S = -1 + 2*(T-T1)/H AND H = T2-T1. C THE COEFFICIENTS ARE... SUBROUTINE GETCO(Y1,Y1P,Y2,Y2P) Implicit DoublePrecision (A-H,O-Z) COMMON /CFTS/ COEF(4),H8 COEF(3)=H8*(Y2P-Y1P) COEF(1)=(Y2+Y1)/2.0 - COEF(3) COEF(4)=H8*(Y2P+Y1P)-(Y2-Y1)/4.0 COEF(2)=(Y2-Y1)/2.0 - COEF(4) RETURN END C C EVALUATE P(S) DoublePrecision FUNCTION EVAL(S) Implicit DoublePrecision (A-H,O-Z) COMMON /CFTS/ COEF(4),H8 EVAL=COEF(1)+S*(COEF(2)+S*(COEF(3)+S*COEF(4))) RETURN END