C #DFFSB V1A 2-MAY-73. C LAST UPDATE 21-MAY-73. C SUBROUTINE DFFSB( DERIV, Y0, Y1, N, X0, X1, IERR, LOG, OPT) DIMENSION Y0(N), Y1(N), OPT(6) C C *PURPOSE. C THIS IS AN EXTRAPOLATED MIDPOINT INTEGRATION OF A SET C OF "N" DIFFERENTIAL EQUATIONS FROM X0 TO X1 WITH AUTOMATIC C STEP CONTROL. IT CALLS AND IN FACT SUPERVISES "DFFSA" C WHICH IS A SLIGHT MODIFICATION OF "DIFFSYS" BY BULLIRSCH C AND STOER. C C *PARAMETERS: C C DERIV - EXTERNAL NAME OF SUBROUTINE PROVIDING THE RIGHT C HAND SIDES OF THE ODE'S. CALLED AS: C CALL DERIV( Y, YP, N, X) C WHERE Y(N) ARE THE DEPENDENT VARIABLES, X IS THE C INDEPENDENT VARIABLE, AND YP(N) MUST RETURN THE C DERIVATIVES OF THE Y. C Y0 - REAL ARRAY OF LENGTH "N" GIVING THE INITIAL Y C VALUES. C Y1 - REAL ARRAY OF LENGTH "N" RETURNING THE FINAL Y C VALUES. Y0 AND Y1 MAY IN FACT BE THE SAME ARRAY. C N - INTEGER NUMBER OF DEPENDENT VARIABLES Y(I) C X0 - REAL INITIAL VALUE OF INDEPENDENT VARIABLE C X1 - REAL FINAL VALUE OF INDEPENDENT VARIABLE. C IERR - ERROR FLAG RETURNED AS: C =0 NO ERRORS SUCCESSFUL RETURN. C >0 INTEGRATION FAILED X1 NOT REACHED. C <0 INPUT PARAMETERS MADE NO SENSE - ABORTED. C LOG - OUTPUT DEVICE FOR LOGGING MESSAGES FROM DFFSB, SET C TO ZERO IF NO LOG REQUIRED. THE LOG OUTPUT INCLUDES C THE VALUES OF "N", "X0", AND "X1", THE NUMBER OF C MACRO STEPS, THE MAX ESTIMATED ERRORS AND THE RANGE OF C MICRO STEPS USED, AS WELL AS ANY ERROR MESSAGES. C OPT - REAL ARRAY OF LENGTH 6 OR MORE CONTAINING CONTROL C PARAMETERS ETC WHICH ALLOW GREATER CONTROL OVER THE C INTEGRATION. IF THESE VALUES ARE SET ZERO OR -VE "DFFSB" C SUPPLIES DEFAULT VALUES WHICH COVER MOST CASES. C OPT(1) - THE MAXIMUM TOLERABLE RELATIVE ERROR IN ANY Y(I) C DEFAULT IS 1.E-6 C OPT(2) - THE MAXIMUM TOLERABLE ABSOLUTE ERROR IN ANY Y(I). C DEFAULT IS 1.E-6 C OPT(3) - THE MAXIMUM (MACRO) STEP SIZE ALLOWED. C DEFAULT IS (X1-X0) C OPT(4) - THE MINIMUM (MICRO) STEP ALLOWED. C DEFAULT IS 1.E-7*(X1-X0) C OPT(5) - THE INITIAL (MACRO) STEP SIZE TO USE - DEFAULT C (X1-X0) - RETURNED AS THE LAST (MACRO) STEP SIZE C USED. C OPT(6) - THE LAST X VALUE SUCCESSFULLY REACHED, USEFUL C WHEN IERR RETURNED = +1. C C INPUT PARAMETERS: C DERIV, Y0, N, X0, X1, LOG, OPT(I),I=1,2,3,4,5 C OUTPUT PARAMETERS: C Y1, IERR, OPT(5), OPT(6) C C *METHOD. C THE EXTRAPOLATED MIDPOINT METHOD HAS BEEN DISCUSSED BY C BULLIRSCH AND STOER ( COMM A.C.M ) AND THE ROUTINE C "DFFSA" IS A SLIGHT ADAPTATION OF THEIR ROUTINE. THE BASIC C IDEA IS TO INTEGRATE OVER ONE MACRO STEP BY TAKING FIRST 2 THEN C 4 THEN 8 ETC MICRO-STEPS, USING A SIMPLE INTEGRATION RULE EACH C TIME AND KEEP EXTRAPOLATING ALL THE RESULTS SO FAR DOWN TO C ZERO STEP SIZE. WHEN THE EXTRAPOLATIONS FOR N AND 2*N C STEPS AGREE SUFFICIENTLY WELL THE RESULT IS TAKEN FOR THE C ENDPOINT OF THE MACRO-STEP. C "DFFSB" CALLS "DFFSA" WHICH IS RESPONSIBLE FOR ONE STEP C AND SUPERVISES ITS OPERATION MAKING SURE THAT IT DOES IN C FACT REACH "X1" IF AT ALL POSSIBLE. C THE ACCURACY CRITERION APPLIED TO THE RESULT OF EACH C MACRO-STEP IS THAT THE VALUES OF Y(I) FOUND SHOULD HAVE ESTIMATED C ERRORS SMALLER THAN "RELTOL"*Y(I) + "ABSTOL" WHERE "RELTOL" C AND "ABSTOL" ARE THE RELATIVE AND ABSOLUTE TOLERABLE ERRORS. C A SPECIAL FEATURE OF THIS ROUTINE IS THAT THE CALLER CAN C MONITOR THE INTEGRATION BY PROVIDING A ROUTINE "DECHK" CALLED AS C CALL DECHK( Y, N, X, STOP) C THIS ROUTINE IS CALLED BY "DFFSB" AFTER EVERY SUCCESSFUL MACRO-STEP C WITH THE LATEST X AND Y VALUES. IF THE ROUTINE "DECHK" C SETS THE LOGICAL VARIABLE "STOP" = .TRUE. THEN AN IMMEDIATE C EXIT IS DONE FROM "DFFSB". THIS ALLOWS THE USER TO STOP THE C INTEGRATION WHEN A Y(I) REACHES A CERTAIN VALUE OR A C DERIVATIVE REACHES A CERTAIN VALUE, OR PRINT OUT RESULTS C WHENEVER THEY CHANGE BY A SIGNIFICANT AMOUNT ETC ETC. THE LIBRARY C PROVIDES A DUMMY ROUTINE "DECHK" WHICH DOES NOTHING, ANY ROUTINE C SUPPLIED BY THE USER WILL SUPPLANT THIS ONE. C C *ACCURACY. C THIS METHOD IS ONE OF THE BEST FOR HIGH ACCURACY INTEGRATIONS C OF DIFFERENTIAL EQUATIONS. HOWEVER IT CAN BE QUITE SENSITIVE C TO THE PRECISION OF THE FLOATING POINT SYSTEM AND CAN EVEN C BE WORSE RELATIVE TO ANOTHER ROUTINE (SUCH AS "RUNGE-KUTTA") C ON A MACHINE WITH SIGNIFICANCE OF 2**(-22) OR HIGHER, WHILE IT C WOULD NORMALLY BE MUCH BETTER ON A HIGH PRECISION SYSTEM. THE C CURRENT ROUTINE IS IN SINGLE PRECISION AND DOUBLE PRECISION C SHOULD BE CONSIDERED. C NORMALLY THE PER-STEP ACCURACY ACHIEVED IS SIGNIFICANTLY C BETTER THAN THE TOLERABLE ERROR, OFTEN BY SEVERAL ORDERS OF C MAGNITUDE. BECAUSE OF THIS IT IS ADVISABLE TO EXAMINE THE C REAL ACCURACY OF SOME RESULTS TO SEE IF IN FACT THE TOLERABLE C ERRORS CAN BE INCREASED. C C *RESTRICTIONS. C C *ERROR CONDITIONS. C SEE DETAILS FOR "IERR" ABOVE. C C *NON STANDARD ROUTINES CALLED. C DERIV - AS MENTIONED ABOVE. C DFFSA - BASIC STEP IN EXTRAPOLATED MIDPOINT INTEGRATION. C DECHK - THE MONITORING ROUTINE , SEE ABOVE. C CPRAY - CALL CPRAY(B, A, N) COPIES B(I),I=1,N TO A(I). C C *TYPICAL TIMES. C C *ORIGIN. C BULLIRSCH AND STOER AS MENTIONED ABOVE VIA M.R.MANNING. C C *COMMENTS. C THE INTERFACE WITH THIS ROUTINE HAS BEEN DELIBERATELY MADE C AS SIMILAR AS POSSIBLE TO THAT WITH "RNKTB" SO THE USER CAN C EASILY CHANGE FROM ONE ROUTINE TO ANOTHER. C C #END. C C #START. LOGICAL NOTEND, SMALL, STOP DIMENSION DETAIL(4) DATA PRNG1, PRNG2, EPS, SIG/ 1.2, 0.5, 1.0E-6, 1.0E-7/ C PRNG1 - RATIO BY WHICH DOUBLE STEP MAY BE TEMPORARILY INCREASED C IN ORDER TO HIT THE OUTPUT POINT EXACTLY. C PRNG2 - RATIO BELOW WHICH THE STEP MAY BE REDUCED TEMPORARILY C BUT WILL NOT COUNT TOWARDS SUBSEQUENT STEP REFINEMENT. C EPS - DEFAULT VALUE FOR RELATIVE AND ABSOLUTE TOLERANCES. C SIG - MACHINE SIGNIFICANCE. C IF (LOG .GT. 0) WRITE(LOG, 908) N, X0, X1 908 FORMAT(' DFFSB: INTEGRATING ',I3,' VARIABLES FROM ',G13.6, * ' TO ',G13.6) C IREF= 0 ICSN= 0 ISTP= 0 DELT= X1-X0 ABSDEL= ABS(DELT) EPSDEL= SIG*ABSDEL X=X0 GMIN= ABSDEL GMAX= 0. ABSERR= 0. RELERR= 0. CALL CPRAY( Y0, Y1, N) C C FETCH IN STEERING PARAMETERS OPT(I) C C RELATIVE TOLERABLE ERROR OPT(1) RELTOL= OPT(1) C ABSOLUTE TOLERABLE ERROR OPT(2) ABSTOL= OPT(2) C MAX (MACRO-)STEP SIZE OPT(3) HMAX= OPT(3) C MIN (MICRO-)STEP SIZE OPT(4) HMIN= OPT(4) C START STEP SIZE OPT(5) PH= OPT(5) C LAST X VALUE REACHED OPT(6) - OUTPUT ONLY. C C SET DEFAULTS: IF(RELTOL .LT. 0.) RELTOL= EPS IF(ABSTOL .LT. 0.) ABSTOL= EPS IF( (RELTOL + ABSTOL) .LE. 0.) RELTOL= EPS IF( HMAX .LE. 0. .OR. HMAX .GT. ABSDEL) HMAX=ABSDEL IF( PH .LE. 0. .OR. PH .GT. HMAX) PH= HMAX IF (HMIN .LE. 0.) HMIN= EPSDEL IF( HMIN .GT. HMAX .OR. DELT .EQ. 0.) GO TO 520 C C RECYCLE TO HERE AFTER EVERY SUCCESSFUL STEP C 115 H= SIGN( PH, DELT) TMP= ABS(X1 - X) NOTEND= TMP .GT. PRNG1*PH IF( NOTEND) GO TO 125 C IF LAST STEP TO ENDPOINT AND RATHER SMALL DONT USE TO C COUNT TOWARDS STEP REFINING. SMALL= TMP .LT. PRNG2*PH H= X1 - X C C DO INTEGRATION STEP C 125 OX= X CALL DFFSA(DERIV, Y1, N, X, H, RELTOL, ABSTOL, HMIN, KFLAG, * DETAIL) IF(KFLAG .LT. 0) GO TO 510 ISTP= ISTP+1 G=DETAIL(2) C PICKED UP MICRO-STEP USED NOW UPDATE MIN-MAX LOG. IF( G .GT. GMIN) GO TO 127 C SET NEW MIN GMIN= G JMIN= DETAIL(1) XMIN= OX 127 IF( G .LT. GMAX) GO TO 128 C SET NEW MAX GMAX= G JMAX= DETAIL(1) XMAX= OX C C MAINTAIN PER-STEP ERROR LOG 128 ABSERR= AMAX1( ABSERR, DETAIL(4) ) RELERR= AMAX1( RELERR, DETAIL(3) ) C C CHECK IF ENDPOINT ACTUALLY REACHED. NOTEND= ABS(X1 - X) .GT. EPSDEL C SET NEW STEP SIZE IF APPROPRIATE IF (SMALL) GO TO 130 PH= AMIN1(ABS(H), HMAX) PH= AMAX1(PH, 8.*HMIN) 130 STOP= .FALSE. CALL DECHK( Y1, N, X, STOP) C IF( NOTEND .AND. (.NOT.STOP) ) GO TO 115 C C GOT TO ENDPOINT -- LEAVE NOW. C IF( ( LOG .GT. 0) .AND. STOP) WRITE(LOG, 904) X 904 FORMAT(' DFFSB: TERMINATED BY USERS "DECHK" AT',G13.6) IERR=0 GO TO 610 C C SET OUTPUT PARAMETERS LOG OUT ETC. C C EXIT WHEN STUCK 510 IF( LOG .GT. 0) WRITE(LOG, 902) X 902 FORMAT(' DFFSB: MICRO-STEP UNDERFLOW AT ',G13.6) IERR=1 GO TO 610 C C EXIT - CALL DOES NOT MAKE SENSE. C 520 IF( LOG .GT. 0) WRITE(LOG, 900) 900 FORMAT(' DFFSB: INPUT PARAMETERS PREVENT INTEGRATION') IERR= -1 C C 610 OPT(5)=ABS(H) OPT(6)=X IF (LOG .GT. 0) WRITE(LOG, 906) ISTP, RELERR, ABSERR, GMIN * ,XMIN, JMIN, GMAX, XMAX, JMAX 906 FORMAT(' DFFSB: MADE ',I4,' MACRO-STEPS'/ *' DFFSB: MAX ESTIMATED ERRORS ',G13.6,' RELATIVE ' 2 ,G13.6,' ABSOLUTE'/ *' DFFSB: MICRO-STEPS FROM ',G13.6,' AT ',G13.6,' EXTRAP LEVEL',I3/ *' DFFSB: TO ',G13.6,' AT ',G13.6,' EXTRAP LEVEL',I3) RETURN END