C #DFFSA V1B APR-73. C LAST UPDATE: 21-MAY-73. C SUBROUTINE DFFSA(F, Y, N, X, H0, RELTOL, ABSTOL, GMIN, * KFLAG, DETAIL) DIMENSION Y(N), DETAIL(4) C C *PURPOSE. C THIS ROUTINE DOES A SINGLE STEP OF THE EXTRAPOLATED C MIDPOINT INTEGRATION FOR A SET OF "N" ORDINARY DIFFERENTIAL C EQUATIONS. THE METHOD IS KNOWN AS THE BULLIRSCH STOER C METHOD. C C *PARAMETERS: C C F - DERIVATIVE ROUTINE USED AS: CALL F( Y, DY, N, X) C TO RETURN DERIVATIVES (DY(I), I= 1, N) AT C ( Y(I), I= 1, N) AND X. C Y - REAL ARRAY OF "N" DEPENDENT VARIABLES AT X. C N - INTEGER NUMBER OF DEPENDENT VARIABLES. C X - REAL INDEPENDENT VARIABLE, INPUT AS INITIAL VALUE, RETURNED C AS FINAL VALUE. C H0 - REAL STEP IN X. C RELTOL - RELATIVE ERROR TO BE TOLERATED. C ABSTOL - ABSOLUTE ERROR TO BE TOLERATED. C GMIN - MINIMUM MICRO-STEPSIZE ALLOWED C KFLAG - EXIT INDICATOR FLAG WITH MEANINGS: C +1 SUCCESSFUL STEP OF INPUT H0 TAKEN. C 0 SUCCESSFUL STEP OF LESS THAN INPUT H0 TAKEN C -1 MICRO STEPS UNDERFLOWED GMIN WITHOUT SUCCESS. C DETAIL - ARRAY RETURNING INTERNAL DETAILS OF STEP: C DETAIL(1): EXTRAPOLATION LEVEL USED. C DETAIL(2): MICRO-STEP (ABSOLUTE VALUE) USED. C DETAIL(3): MAX ESTIMATED REL ERROR C DETAIL(4): MAX ESTIMATED ABS ERROR C C INPUT PARAMETERS: C F, Y, N, X, H0, RELTOL, ABSTOL, GMIN C OUTPUT PARAMETERS: C Y, X, H0, KFLAG, DETAIL C C *METHOD. C ONE MACRO-STEP OF THE EXTRAPOLATED MIDPOINT METHOD C IS TAKEN. THE STEP MOVED THROUGH WILL BE EITHER H0 OR H0/2 C OR H0/4 ETC. TO ACHIEVE THE REQUIRED ACCURACY. THE ARRAY "Y" C AND "X" ARE UPDATED TO THE VALUES AT THE END OF THE STEP AND C "H0" IS RETURNED AS THE SUGGESTED STEP SIZE FOR THE NEXT C MACRO-STEP. THE STEPS ARE CHOSEN SO THAT NO MORE THAN 9 C EXTRAPOLATION LEVELS SHOULD BE NECESSARY FOR THE NEXT STEP. C THE ACCURACY CRITERION IS THAT ALL THE COMPONENTS Y(I) C SHOULD HAVE ESTIMATED ERRORS SMALLER THAN "RELTOL"*Y(I) + "ABSTOL" C C C *ACCURACY. C SEE COMMENTS TO "DFFSB" C C *RESTRICTIONS. C CURRENT VERSION LIMITED TO N<=20 BY INTERNAL DIMENSIONS.! C C *ERROR CONDITIONS. C SEE DETAILS OF "KFLAG" ABOVE. C C *NON STANDARD ROUTINES CALLED. C F USER SUPPLIED DERIVATIVE ROUTINE. C C *TYPICAL TIMES. C C *ORIGIN. C R.BULLIRSCH AND J.STOER NUMERISCHE MATHEMATIK 8 (66) 1. C TRANSLATED INTO FORTRAN AND THE ACCURACY CRITERION CHANGED, C AS WELL AS GIVING THE CALLING ROUTINE MORE CONTROL AND INFORMATION. C C *COMMENTS. C C #END. C C #START. C INTEGER R, SR DIMENSION YA(20), YL(20), YM(20), DY(20), DZ(20), DT(20,7) DIMENSION D(7), YG(8,20), YH(8,20) LOGICAL KONV, BO, BH CSTART CONA=16./9. CONB= 64./9. CONC= 256./9. CALL F( Y, DZ, N, X) BH= .FALSE. KFLAG= 1 DO 10 I= 1, N 10 YA(I)= Y(I) C CA-RECYCLE TO HERE: 50 A= H0 + X FC= 1.5 BO= .FALSE. M= 1 R= 2 SR= 3 JJ= 0 DO 500 J= 1, 10 CLOOP-EXTRAPOLATION LEVEL IF (BO) GO TO 60 CIFNOT-BO D(2)= 2.25 D(4)= 9. D(6)= 36. GO TO 70 CELSIF-BO 60 D(2)= CONA D(4)= CONB D(6)= CONC CENDIF-BO C 70 KONV= J .GT. 3 IF ( J .GT. 7) GO TO 80 CIFNOT-J > 7 L= J D(L)= M*M GO TO 90 CELSIF-J > 7 80 L= 7 D(7)= 64. FC= 0.6*FC CENDIF-J C 90 M= M + M G= H0/FLOAT(M) IF (G .LT. GMIN) GO TO 650 B= G + G IF ( BH .AND. ( J .LT. 9) ) GO TO 200 CIFNOT-MAJOR KK= ( M - 2)/2 M= M - 1 DO 110 I= 1, N YL(I)= YA(I) 110 YM(I)= YA(I) + G*DZ(I) DO 150 K= 1, M CALL F( YM, DY, N, X + FLOAT(K)*G) DO 120 I= 1, N U= YL(I) + B*DY(I) YL(I)= YM(I) 120 YM(I)= U IF ( ( K .NE. KK) .OR. ( K .EQ. 2) ) GO TO 150 CIFNOT-K JJ= JJ + 1 DO 130 I= 1, N YH(JJ, I)= YM(I) 130 YG(JJ, I)= YL(I) CENDIF-K 150 CONTINUE GO TO 300 CELSIF-MAJOR 200 DO 210 I= 1, N YM(I)= YH(J, I) 210 YL(I)= YG(J, I) CENDIF-MAJOR 300 CALL F( YM, DY, N, A) C CLEAR ESTIMATES OF ERRORS FOR LOGGING PURPOSES ABSERR= 0. RELERR= 0. DO 400 I= 1, N V= DT(I, 1) TA= ( YM(I) + YL(I) + G*DY(I) )/2. C= TA DT(I, 1)= TA DO 350 K= 2, L B1= D(K)*V B= B1 - C U= V IF ( B .EQ. 0.) GO TO 320 CIFNOT-B B= (C - V)/B U= C*B C= B1*B CENDIF-B 320 V= DT(I, K) DT(I, K)= U 350 TA= TA + U U= ABS(TA) V= ABS(Y(I) - TA) C UPDATE MAX PER-STEP ERROR ESTIMATES ABSERR= AMAX1( ABSERR, V) IF (TA .NE. 0.) RELERR= AMAX1(RELERR, V/U ) IF (V .GT. (RELTOL*U + ABSTOL) ) KONV= .FALSE. 400 Y(I)= TA IF (KONV) GO TO 700 D(3)= 4. D(5)= 16. BO= .NOT. BO M= R R= SR 500 SR= M + M CENDO-EXTRAPOLATION LEVEL BH= .NOT. BH KFLAG= 0 H0= H0/2. GO TO 50 CA-RECYCLE C C SET ERROR RETURN CODE 650 KFLAG= -1 C C EXIT ON CONVERGENCE. 700 H0= FC*H0 C RETURN EXTRAPOLATION LEVEL DETAIL(1)= J C RETURN MICRO STEP DETAIL(2)= ABS(G) C RETURN LAST ESTIMATED MAX PER-STEP ERRORS DETAIL(3)= RELERR DETAIL(4)= ABSERR X= A RETURN END