C #RNKTB V1C 19-JAN-73. C LAST UPDATE: 30-MAY-73. C SUBROUTINE RNKTB(DERIV,Y0,Y1,NN,X0,X1,IERR,LOG,OPT,WK) DIMENSION Y0(NN), Y1(NN), OPT(10), WK(NN,9) C C *PURPOSE. C TO INTEGRATE A SET OF "NN" FIRST ORDER ORDINARY DIFFERENTIAL C EQUATIONS USING A VARIABLE STEP RUNGE-KUTTA FOURTH-ORDER METHOD. C THE INTEGRATION IS FROM X0, Y0(I) TO X1, Y1(I). C C *PARAMETERS. C C DERIV- IS THE EXTERNAL NAME OF THE USER'S SUBROUTINE GIVING THE C RIGHT HAND SIDES OF THE DIFFERENTIAL EQUATIONS. IT IS C CALLED IN THE FORM: C CALL DERIV(Y,YP,NN,X) C WHERE Y GIVES THE VALUES OF THE "NN" DEPENDENT VARIABLES C AND X IS THE INDEPENDENT VARIABLE. THE ROUTINE DERIV C MUST RETURN THE DERIVATIVES OF THE Y(I) IN THE ARRAY C YP OF LENGTH NN. C Y0 - IS A REAL ARRAY OF LENGTH NN GIVING THE INITIAL Y VALUES C Y1 - IS A REAL ARRAY OF LENGTH NN RETURNING THE FINAL Y VALUES C Y0 AND Y1 MAY IN FACT BE THE SAME USER'S ARRAY. C NN - IS THE NUMBER OF EQUATIONS BEING INTEGRATED. C X0 - REAL VALUE, INTIAL VALUE OF INDEPENDENT VARIABLE X. C X1 - REAL VALUE, FINAL VALUE OF INDEPENDENT VARIABLE X. C IERR - IS AN ERROR FLAG. RETURNED ZERO IF THE INTEGRATION SUCCESSFULLY C REACHED X1. IS (-1) IF AN INPUT PARAMETER WAS FOUND C TO BE IMPROPER. AND IS +VE IF THE INTEGRATION FAILED TO C REACH X1 (AND OPT(7) IS -VE) IN WHICH CASE IERR C IS THE INDEX OF THE FIRST Y COMPONENT FOR WHICH THE ERROR C COULD NOT BE MADE SMALL ENOUGH. C LOG - IS THE DEVICE NUMBER FOR LOGGING INFORMATION FROM "RNKTB" C OR IS ZERO IF NO PRINTOUT REQUIRED FROM "RNKTB". THE C LOG OUTPUT IS CLEARLY IDENTIFIED AS COMING FROM "RNKTB" C AND GIVES VALUES OF "NN", "X0", "X1" TOGETHER WITH C NUMBER OF DERIVATIVE CALCULATIONS, STEPS TAKEN, STEP SIZE C CHANGES AND MIN AND MAX STEP SIZES. C IN THE EVENT OF AN ERROR A MESSAGE IS ALSO OUTPUT. C OPT - IS A REAL ARRAY OF DIMENSION 10 OR MORE CONTAINING CONTROL C PARAMETERS ETC. WHICH ALLOW GREATER CONTROL OVER C "RNKTB". WHILE NOT NECESSARY TO SET ANY OF THESE IT MAY C HELP IN SOME CIRCUMSTANCES. THE DEFAULT VALUES OF C INPUT "OPT" VALUES ARE INVOKED BY LEAVING THEM ZERO. C OPT(1) - IS THE MAX ALLOWED RELATIVE ERROR - DEFAULT 1.E-6 C OPT(2) - IS THE MAX ALLOWED ABSOLUTE ERROR - DEFAULT 1.E-6 C OPT(3) - IS THE MAX ALLOWED STEP SIZE - DEFAULT (X1-X0)/2 C OPT(4) - THE MIN ALLOWED STEP SIZE - DEFAULT (X1-X0)/2**16 C OPT(5) - IS SET INITIALLY TO THE INITIAL STEP SIZE TO BE C USED - DEFAULT IS (X1-X0)/64 - AND RETURNED AS C LAST STEP SIZE USED. C OPT(6) - IS SET BY "RNKTB" TO BE THE LAST X VALUE REACHED C SUCCESSFULLY, USEFUL WHEN IERR IS RETURNED =1. C OPT(7) - IS SET +VE TO FORCE "RNKTB" TO KEEP GOING WITH C MINIMUM SIZED STEPS EVEN WHEN THE TOLERABLE ERROR C CANNOT BE ACHIEVED - DEFAULT IS ZERO OR NEGATIVE C CAUSING TERMINATION OF THE INTEGRATION. C OPT(8) - IS THE FACTOR BETWEEN THE MAXIMUM TOLERABLE ERROR C AND THE ERROR AT WHICH STEP SIZE INCREASE SHOULD C BE CONSIDERED - DEFAULT IS 32. C OPT(9) - IS THE NUMBER OF STEPS WITH "TOO-GOOD" ERRORS C WHICH MUST BE TAKEN BEFORE THE STEP SIZE IS INCREASED C - DEFAULT IS 1. C OPT(10)- IS RETURNED AS THE NUMBER OF STEPS WHICH HAVE C JUST BEEN TAKEN AND WERE "TOO-GOOD". C WK - IS AN ARRAY OF AT LEAST 9*NN ELEMENTS USED AS WORK SPACE C BY THE ROUTINE. C C INPUT PARAMETERS: C DERIV, Y0, NN, X0, X1, LOG, OPT(I):I=1 TO 5, AND 7 TO 10. C OUTPUT PARAMETERS: C Y1, IERR, OPT(I):I=5, 6, AND 10, WK C C C *METHOD. C THE METHOD IS BASED ON THE FOURTH ORDER RUNGE-KUTTA ALGORITHM C TAKEN OVER ONE STEP SIZE AND ANOTHER STEP DOUBLE THE SIZE. THE C DIFFERENCE IN FINAL VALUES PREDICTED FOR THE TWO STEP SIZES GIVES C A MEASURE OF THE ERROR AND ALLOWS THE ROUTINE TO REDUCE OR INCREASE C THE STEP SIZE IF NECESSARY. THE RESULTS FOR THE TWO STEP SIZES C ALSO ALLOW A ROUGH EXTRAPOLATION TO A RESULT FOR ZERO STEP SIZE C (RICHARDSON EXTRAPOLATION) AND PROVIDED THE STEP SIZES ARE NOT TOO C LARGE THIS GENERALLY IMPROVES THE ACCURACY BY AN ORDER OF MAGNITUDE C AND MAKES THE ACTUAL ACCURACY BETTER THAN THE TOLERABLE ERRORS C IN OPT(1) AND OPT(2). C THE ACCURACY CRITERION APPLIED TO EACH STEP IS THAT THE C RESULT SHOULD HAVE AN ESTIMATED ERROR OF LESS THAN C "RELTOL"*Y(I) + "ABSTOL" IN EACH Y(I). WHERE RELTOL AND C ABSTOL ARE THE MAXIMUM TOLERABLE RELATIVE AND ABSOLUTE C ERRORS. C AS A SPECIAL FEATURE OF THIS ROUTINE THE CALLER CAN MONITOR C THE INTEGRATION IN DETAIL BY PROVIDING A ROUTINE "DECHK" CALLED AS: C CALL DECHK(Y,NN,X,STOP) C THIS ROUTINE IS CALLED AFTER EVERY SUCCESSFUL STEP BY "RNKTB" WITH C THE LATEST X AND Y VALUES. IF THIS ROUTINE SETS LOGICAL STOP=.TRUE. C "RNKTB" WILL RETURN IMMEDIATELY TO ITS CALLING PROGRAM. THIS ALLOWS C THE USER TO INTEGRATE UNTIL A Y(I) VALUE REACHES SOME CRITICAL C VALUE. OR UNTIL CHANGES IN Y(I) OR ITS DERIVATIVE ARE LESS THAN C SOME CRITICAL VALUE ETC. ETC. THE LIBRARY PROVIDES A DUMMY "DECHK" C ROUTINE WHICH DOES NOTHING, ANY ROUTINE SUPPLIED BY THE USER WITH C THIS NAME WILL SUPPLANT THE DUMMY ONE. C C *ACCURACY. C THE ROUTINE ATTEMPTS TO INTEGRATE TO A SPECIFIED ACCURACY. C IN SOME CASES OF UNSTABLE EQUATIONS THE ROUTINE WILL C THINK IT HAS ACHEIVED THE REQUIRED ACCURACY WHEN IT HAS NOT. C HOWEVER, GENERALLY THE ROUTINE ACHIEVES A BETTER PER-STEP ACCURACY C THAN SPECIFIED, OR STOPS, OR IS FORCED TO CONTINUE BY OPT(7) +VE. C C *RESTRICTIONS. C C *ERROR CONDITIONS. C IF A LOGGING DEVICE IS SPECIFIED "RNKTB" WILL GIVE ADEQUATE ERROR C MESSAGES. IF NO LOGGING DEVICE IS SPECIFIED THEN THE CALLER SHOULD C CHECK "IERR" ON RETURN AND IN MANY CASES IT WOULD BE AS WELL TO C CHECK "IERR" ON RETURN ANYWAY. THE MEANINGS OF NON-ZERO IERR ARE C GIVEN ABOVE. C C *NON-STANDARD ROUTINES CALLED. C DERIV(USER'S) C CPRAY, RKSTB, (MRMLIB) C DECHK(USER OR MRMLIB). C C *TYPICAL TIMES. C C *ORIGIN. C "RNKTB" OWES MUCH TO A SIMILAR ROUTINE "RUNKUT" BY R.E.JONES C SANDIA LABORATORIES, ALBUQUERQUE. THAT ROUTINE IS THOROUGHLY C DOCUMENTED WITH EXAMPLES IN THE REPORT SC-M-70-724, NOV-1970. C "RNKTB" DOES NOT OFFER SOME OF THE FEATURES OF "RUNKUT" - NOTABLY C THE STEP SIZES ARE NOW EITHER DOUBLED OR HALVED. C C *COMMENTS. C VARIABLE STEP SIZE INTEGRATORS FOR DIFFERENTIAL EQUATIONS C ARE CURRENTLY THE SUBJECT OF MUCH INVESTIGATION. THE ROUTINE C "RNKTB" IS BY NO MEANS THE MOST EFFICIENT KNOWN BUT IT IS A VERY C RELIABLE ROUTINE. THE METHOD OF BULLIRSCH AND STOER IS PROBABLY C A BETTER ROUTINE IN GENERAL AND CERTAINLY ANY ROUTINE WHICH IS C TO INTEGRATE OVER LARGE INTERVALS WITH RELATIVELY FEW STEP SIZE C CHANGES MUST DO BETTER THAN HALVING OR DOUBLING THE STEP! C C #END. C C LOGICAL NOTEND, SMALL, STOP C INDICES RELATING WK COLUMNS TO VECTORS ARE: DATA NYP,NYY,NYPA,NDY1,NDYA,NDY2,NER,NYA/1,2,3,4,5,6,7,9/ C DATA PRNG1, PRNG2, EPS/ 2.1, 1.9, 1.0E-6/ C PRNG1- RATIO BY WHICH DOUBLE STEP MAY BE TEMPORARILY INCREASED C IN ORDER TO HIT THE ENDPOINT EXACTLY. C PRNG2- RATIO BELOW WHICH STEP TO OUTPUT POINT IS NOT COUNTED C TOWARDS STEP REFINING. C EPS - DEFAULT VALUE FOR RELATIVE AND ABSOLUTE TOLERANCES. C C #START HERE:: INITIALIZE EVERYTHING IN SIGHT: C IF(LOG .GT. 0) WRITE(LOG, 908) NN, X0, X1 908 FORMAT(' RNKTB: INTEGRATING ',I3,' VARIABLES FROM ',G13.6, * ' TO ',G13.6) C IDCAL=0 IREF=0 ICSN=0 ISUC=0 DELT= X1 - X0 X= X0 SMALL= .FALSE. PHMIN= ABS(DELT) PHMAX= 0. CALL CPRAY( Y0, Y1, NN) C C FETCH STEERING PARAMETERS OPT(I) C C MAX RELATIVE ERROR RELMAX=OPT(1) C MAX ABSOLUTE ERROR ABSMAX=OPT(2) C MAX STEP SIZE HMAX=OPT(3) C MIN STEP SIZE HMIN=OPT(4) C STARTING STEP SIZE -- OUTPUT IS CURRENT STEP SIZE! PH=OPT(5) C LAST X VALUE SUCCESSFULLY REACHED - OPT(6) - OUTPUT ONLY. C STEP CONTROL INDICATOR -- OPT(7) C LE.0. MEANS FAIL IF STEPS GO TOO SMALL ELSE MAINTAIN MINIMUM C STEP SIZE WHEN .GE.0. C FACTOR BETWEEN ALLOWED ERROR AND TOO-GOOD ERROR OPT(8) FAC= OPT(8) C NUMBER OF T00-GOOD STEPS BEFORE INCREASING STEP SIZE NTG= OPT(9) C NUMBER OF JUST PAST TOO-GOOD STEPS ITG= OPT(10) C C SET DEFAULTS AS NECESSARY. IF(RELMAX.LT.0.) RELMAX=EPS IF(ABSMAX.LT.0.) ABSMAX=EPS IF(RELMAX+ABSMAX.LE.0.) RELMAX= EPS TMP=ABS(DELT)/2. IF(HMAX.LE.0.) HMAX=TMP IF(PH.LE.0.) PH=AMAX1(HMAX/32.,HMIN) IF(PH.GT.HMAX) PH=HMAX IF(HMIN.LE.0.) HMIN=PH/1024. IF(HMIN.GT.HMAX .OR. DELT.EQ.0.) GO TO 520 IF(FAC .LE. 0.) FAC= 32. IF(NTG .LE. 0) NTG= 1 IF(ITG .LT. 0) ITG= 0 C C RECYCLE TO HERE AFTER EVERY SUCCESSFUL INTEGRATION STEP. C C CHECK IF WE CAN REACH END POINT X1. C 115 H=SIGN(PH,DELT) TMP=ABS(X1-X) NOTEND= TMP .GT. PRNG1*PH IF( NOTEND ) GO TO 125 C IF WE ARE DOING LAST STEP SEE IF IT IS TOO SMALL TO COUNT C TOWARDS REFINING RUNNING STEP SIZE. SMALL= TMP .LT. PRNG2*PH H=0.5*(X1-X) C C************************************* C START INTEGRATION STEP. C C FIRST DO DOUBLE-SIZE STEP. C 125 XA=X CALL RKSTB(DERIV,Y1,WK(1,NYP),WK(1,NYY),NN,XA,H+H,0 * ,WK(1,NDY1)) IDCAL=IDCAL+4 C C RECYCLE TO HERE AFTER REDUCING STEP SIZE. C C NOW DO TWO STEPS OF SIZE H C 155 XA=X CALL RKSTB(DERIV,Y1,WK(1,NYP),WK(1,NYA),NN,XA,H,1 * ,WK(1,NDYA) ) C CALL RKSTB(DERIV,WK(1,NYA),WK(1,NYPA),WK(1,NYY),NN,XA,H,0 * ,WK(1,NDY2) ) IDCAL=IDCAL+7 C C************************************* C STEP SIZE CONTROL AREA C C ESTIMATE ERRORS AND COMPARE WITH ALLOWED ERRORS C C REMEMBER PAST T00-GOOD STEPS + CLEAR ITG ITMP= ITG ITG= 0 DO 205 I=1,NN WK(I,NDY2)=WK(I,NDY2) + WK(I,NDYA) C STORE ESTIMATED ERRORS WK(I,NYPA)=ABS(WK(I,NDY1)-WK(I,NDY2) )/30. C STORE ALLOWED ERRORS. WK(I,NER)=RELMAX*ABS(WK(I,NYY) ) + ABSMAX IERR=I IF(WK(I,NYPA).GT.WK(I,NER) ) GO TO 225 205 CONTINUE C DROP THROUGH IF STEP SUCCESSFUL ISUC=ISUC+1 C UPDATE MIN-MAX STEP RECORD. IF( PH.GT.PHMIN) GO TO 2052 C SET NEW MIN STEP -- PHMIN=PH XMIN= X 2052 IF( PH .LT. PHMAX) GO TO 2054 C SET NEW MAX STEP PHMAX= PH XMAX= X 2054 IF( PH .GE. HMAX .OR. SMALL) GO TO 250 C C CHECK IF ALL ERRORS "TOO SMALL" C DO 210 I=1,NN IF(WK(I,NYPA).GT.WK(I,NER)/FAC ) GO TO 250 210 CONTINUE C DROP THROUGH WHEN THIS POINT IS "TOO GOOD". C C COARSEN STEP. C ICSN=ICSN+1 C RECOVER PAST TOO-GOOD STEPS AND INCREMENT ITG= ITMP + 1 IF(ITG .LT. NTG) GO TO 250 C OK HAVE HAD REQUIRED NR OF TOO-GOOD STEPS: C RESET COUNTER ITG AND DOUBLE STEP SIZE. ITG= 0 PH=PH+PH IF(PH.GT.HMAX) PH=HMAX GO TO 250 C C REFINE STEP. C 225 IREF=IREF+1 PH=0.5*ABS(H) IF(PH-HMIN) 226,226,227 C IF STEP SIZE TOO SMALL CHECK OPT(7) FOR INSTRUCTIONS. 226 IF(OPT(7)) 510,510,250 C 227 H=SIGN(PH,DELT) 228 CALL CPRAY(WK(1,NDYA), WK(1,NDY1), NN) C GO TO 155 C C***************************************** C C IF ERRORS ARE O.K., FINISH THE BOOKEEPING, ETC. FOR THIS STEP. C C DO RICHARDSON EXTRAPOLATION TO ZERO STEP SIZE. 250 DO 270 I=1,NN 270 Y1(I)=WK(I,NYY) + (WK(I,NDY2)-WK(I,NDY1) )/15. X=XA C C CHECK FOR OUTPUT CONDITIONS, ETC. C 275 STOP= .FALSE. CALL DECHK( Y1, NN, X, STOP) C 350 IF(NOTEND .AND. (.NOT. STOP) ) GO TO 115 C NOW GOT TO REQUIRED POINT -- OUTPUT RESULTS. IERR=0 IF( (LOG.GT.0) .AND. STOP ) WRITE(LOG,904) X 904 FORMAT(' RNKTB: TERMINATED BY USERS "DECHK" AT',G13.6) GO TO 610 C C SET OUTPUT CODES ETC., WRITE MSG IF NEEDED AND EXIT. C C EXIT--BECAUSE WE ARE STUCK! 510 IF(LOG.GT.0) WRITE(LOG,902) X,IERR 902 FORMAT(' RNKTB: STUCK AT ',G13.6,' ON INDEPENDENT VARIABLE' * ,I3) GO TO 610 C C EXIT--BECAUSE CALL DOESN'T MAKE SENSE! 520 IERR=-1 IF(LOG.GT.0) WRITE(LOG,900) 900 FORMAT(' RNKTB: INPUT PARAMETERS PREVENT INTEGRATION') C C EXIT -- GENERAL. 610 OPT(5)=PH OPT(6)=X OPT(10)= ITG IF(LOG.GT.0) WRITE(LOG,906) IDCAL,ISUC,IREF,ICSN,PHMIN,XMIN * ,PHMAX,XMAX 906 FORMAT(' RNKTB: MADE ',I4,' DERIVATIVE CALLS FOR ',I4, *' SUCCESSFUL DOUBLE STEPS'/ *' RNKTB: MADE ',I4,' STEP REDUCTIONS AND ',I4,' INCREASES'/ *' RNKTB: USED STEP SIZES FROM ',G13.6,' AT ',G13.6/ * ' RNKTB:',19X,'TO ',G13.6,' AT ',G13.6/) RETURN C END