C ALG. 669R, REMARK ON ALG. 669, UNOFFICIAL RELEASE BY RENKA, 5/25/92 C ALGORITHM 669, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 15, NO. 1, PP. 29-30. C PROGRAM RUN (INPUT,OUTPUT,RESULT,TAPE5=INPUT,TAPE6=RESULT) C C PROGRAM TO INTEGRATE DY/DT=-Y, Y(0)=1 IN THE RANGE 0.LE.T.LE.20 C WITH OUTPUT REQUIRED AT T=1.0,2.0,.....,19.0,20.0. C C THE DIMENSION STATEMENT HAS BEEN SET UP FOR 10 DIFFERENTIAL C EQUATIONS. C DIMENSION Y(10),YP(10),WORK(96),IWORK(5) COMMON/VALS/PI,GAM EXTERNAL FCN C PI = 4.0E0*ATAN(1.0E0) GAM = 1.4 C C SET INPUT PARAMETERS TO SUBROUTINE BRKF45 C NEQN = 1 Y(1) = 1.0E+0 T = 0.0E+0 IFLAG = -1 DEL = 1.0E+0 K = 1 TOUT = FLOAT(K)*DEL TEND = 20.0E+0 RELERR = 1.E-10 ABSERR = 1.E-6 WRITE(6, 50) C C INPUT PARAMETERS SET. C 10 CALL BRKF45 (FCN,NEQN,Y,T,TOUT,TEND,YP,RELERR,ABSERR,IFLAG,WORK, * IWORK) IF (IFLAG.LE.2) GO TO 30 C C IF IFLAG IS GREATER THAN 2 WE HAVE AN INPUT ERROR. C WRITE (6,20) IFLAG,NEQN,IWORK(1),T,TOUT,RELERR,ABSERR 20 FORMAT (1X,'FAILURE',5X,'IFLAG,NEQN,NUM F EVAL = ',2I5,I10/ * 13X,'T,TOUT = ',2F13.4,/13X,'RELERR,ABSERR = ',1P2E10.1) IF (IFLAG .EQ. 3) THEN WRITE (6,25) 25 FORMAT (1X,'RELATIVE ERROR TOLERANCE INCREASED TO CONTINUE') GO TO 10 ELSE STOP END IF 30 CONTINUE IF (T.LT.TOUT) GO TO 10 C C THE INTEGRATION HAS PASSED THE OUTPUT POINT T=TOUT. C A = EXP(-T) B = -A E1 = Y(1)-A E2 = YP(1)-B C C E1 AND E2 ARE THE ERRORS IN THE INTERPOLATED SOLUTION C AND DERIVATIVE C WRITE (6,60) T,Y(1),YP(1),E1,E2 C C SET NEW VALUE OF THE OUTPUT POINT. C K = K + 1 TOUT = FLOAT(K)*DEL IF(ABS(TOUT-TEND).LT.ABSERR) TOUT=TEND IF (T.NE.TEND) GO TO 10 50 FORMAT(5X,1HT, 11X, 6H Y, 14X, * 10H DY/DT, 6X, 10HERROR IN Y, 3X, * 14HERROR IN DY/DT) 60 FORMAT (1X,F8.4,1P2E22.10,2E13.2) END SUBROUTINE FCN (T,Y,YP) DIMENSION Y(1),YP(1) YP(1) = -Y(1) RETURN END SUBROUTINE BRKF45 (FCN,NEQN,Y,T,TOUT,TEND,YP,RELERR,ABSERR,IFLAG, * WORK,IWORK) C C BLOCK FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD ADVANCING C A BLOCK OF TWO EQUAL STEPS. THE FORMULA AT THE FIRST STEP IS C 5(4) WHILE AT THE SECOND STEP IT IS 6(4). C C WRITTEN BY J.R.CASH, C DEPARTMENT OF MATHEMATICS, C IMPERIAL COLLEGE, C SOUTH KENSINGTON, LONDON SW7 2AZ, C ENGLAND. C C THIS IS A HEAVILY REVISED VERSION OF RKF45 OF L.F. SHAMPINE AND C H.A. WATTS. C BRKF45 IS PRIMARILY DESIGNED TO SOLVE NON-STIFF AND MILDLY STIFF C INITIAL VALUE ORDINARY DIFFERENTIAL EQUATIONS WHEN DERIVATIVE C EVALUATIONS ARE INEXPENSIVE. BRKF45 USES INTERPOLATION TO PRODUCE C OUTPUT AT 'OFF-STEP POINTS' EFFICIENTLY. BRKF45 SHOULD GENERALLY C NOT BE USED WHEN THE USER IS DEMANDING HIGH ACCURACY. IN SUCH C CASES A GOOD ADAMS CODE WILL OFTEN BE MORE EFFICIENT. C C*********************************************************************** C ABSTRACT C*********************************************************************** C C SUBROUTINE BRKF45 INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C C DY(I)/DT = FCN(T,Y(1),Y(2),...,Y(NEQN)) C C WHERE THE Y(I) ARE KNOWN AT TIME T. C TYPICALLY THE SUBROUTINE IS USED TO INTEGRATE FROM T TO TEND C (WHILE RETURNING ANSWERS AT SPECIFIED OUTPUT POINTS TOUT), BUT C IT CAN ALSO BE USED AS A ONE-BLOCK INTEGRATOR TO ADVANCE THE C SOLUTION A SINGLE BLOCK STEP IN THE DIRECTION OF TEND. ON RETURN C THE PARAMETERS IN THE CALL LIST ARE SET FOR CONTINUING THE C INTEGRATION. THE USER HAS ONLY TO CALL BRKF45 AGAIN (AND PERHAPS C DEFINE A NEW VALUE FOR TOUT). ACTUALLY, BRKF45 IS AN INTERFACING C ROUTINE WHICH CALLS SUBROUTINE RKFC FOR THE SOLUTION. SUBROUTINE C RKFC COMPUTES AN APPROXIMATE SOLUTION OVER ONE BLOCK OF LENGTH 2H. C BRKF45 IS PARTICULARLY USEFUL WHEN OUTPUT IS REQUIRED AT MANY C 'OFF-STEP POINTS' SINCE THE OUTPUT VALUES CAN BE OBTAINED BY C INTERPOLATION. THIS IS IN CONTRAST TO MANY OTHER RUNGE-KUTTA C PROGRAMS WHICH CHOOSE THE STEP SEQUENCE SO AS TO HIT ALL OUTPUT C POINTS EXACTLY AND SO BECOME INEFFICIENT WHEN OUTPUT IS REQUIRED C AT MANY POINTS WITHIN A STEP. C BRKF45 USES THE (5,4); (6,4) BLOCK FORMULA DESCRIBED IN C J.R. CASH, A BLOCK 6(4) RUNGE-KUTTA FORMULA FOR NON-STIFF C INITIAL VALUE PROBLEMS, C TO APPEAR IN TOMS. C C C C THE PARAMETERS REPRESENT- C FCN -- SUBROUTINE FCN(T,Y,YP) TO EVALUATE DERIVATIVES C YP(I)=DY(I)/DT C NEQN -- NUMBER OF DIFFERENTIAL EQUATIONS TO BE INTEGRATED. C Y(*) -- APPROXIMATION TO THE SOLUTION VECTOR AT T. C T -- INDEPENDENT VARIABLE. C TOUT -- THE NEXT POINT WHERE INTERMEDIATE OUTPUT IS C REQUIRED. APPROXIMATE SOLUTIONS AT THESE POINTS C WILL BE OBAINED BY INTERPOLATION. C TEND -- END OF THE INTEGRATION RANGE. THIS WILL BE HIT EXACTLY. C YP(*) -- APPROXIMATION TO THE DERIVATIVE VECTOR DY/DT AT T. C RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR C LOCAL ERROR TEST. AT THE NTH POINT OF EACH BLOCK (N=1,2) C THE CODE REQUIRES THAT C ABS(LOCAL ERROR)/N .LE. RELERR*ABS(Y) + ABSERR C FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION VECTORS C IFLAG -- INDICATOR FOR STATUS OF INTEGRATION. C WORK(*) -- ARRAY TO HOLD INFORMATION INTERNAL TO BRKF45 WHICH C IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE DIMENSIONED C AT LEAST 6+9*NEQN C IWORK(*) -- INTEGER ARRAY USED TO HOLD INFORMATION INTERNAL TO C BRKF45 WHICH IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE C DIMENSIONED AT LEAST 5 C C C*********************************************************************** C FIRST CALL TO BRKF45 C*********************************************************************** C C SUBROUTINE RKFC CONTAINS A MACHINE-DEPENDENT PARAMETER U. C THE SINGLE-PRECISION TOMS VERSION OF THIS PACKAGE CONTAINS A VALUE C SUITABLE FOR A CDC 6600, AND THE DOUBLE-PRECISION VERSION CONTAINS C A VALUE SUITABLE FOR AN IBM 3081. IF THE PACKAGE IS RUN ON ANY C OTHER MACHINE, THE USER SHOULD FIRST RESET U, WHICH IS DEFINED IN C A DATA STATEMENT, CLEARLY MARKED IN SUBROUTINE RKFC. APPROPRIATE C VALUES OF U FOR OTHER MAJOR MACHINES ARE SPECIFIED IN THE COMMENTS C OF RKFC. C C THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE ARRAYS C IN THE CALL LIST - Y(NEQN) , YP(NEQN) , WORK(6+9*NEQN) C IWORK(5), DECLARE FCN IN AN EXTERNAL STATEMENT, SUPPLY SUBROUTINE C FCN(T,Y,YP) AND INITIALIZE THE FOLLOWING PARAMETERS- C C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED. (NEQN .GE. 1) C Y(*) -- VECTOR OF INITIAL CONDITIONS. C T -- STARTING POINT OF INTEGRATION , MUST BE A VARIABLE. C TOUT -- OUTPUT POINT AT WHICH SOLUTION, AND POSSIBLY THE C DERIVATIVE, IS DESIRED. C TEND -- END OF THE RANGE OF INTEGRATION. IF THE SOLUTION IS C REQUIRED ONLY AT TEND THEN THE USER SHOULD SET TOUT=TEND. C RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES C WHICH MUST BE NON-NEGATIVE. RELERR MUST BE A VARIABLE WHILE C ABSERR MAY BE A CONSTANT. THE CODE SHOULD NORMALLY NOT BE C USED WITH RELATIVE ERROR CONTROL SMALLER THAN ABOUT 1.E-8, C UNLESS AN APPROPRIATE NONZERO ABSOLUTE TOLERANCE IS GIVEN. C TO AVOID LIMITING PRECISION DIFFICULTIES THE CODE REQUIRES C RELERR TO BE LARGER THAN AN INTERNALLY COMPUTED RELATIVE C ERROR PARAMETER WHICH IS MACHINE DEPENDENT. IN PARTICULAR, C PURE ABSOLUTE ERROR IS NOT PERMITTED. IF A SMALLER THAN C ALLOWABLE VALUE OF RELERR IS ATTEMPTED, BRKF45 INCREASES C RELERR APPROPRIATELY AND RETURNS CONTROL TO THE USER BEFORE C CONTINUING THE INTEGRATION. C IFLAG -- +1,-1 INDICATOR TO INITIALIZE THE CODE FOR EACH NEW C PROBLEM. NORMAL INPUT IS +1. THE USER SHOULD SET IFLAG=-1 C ONLY WHEN ONE-BLOCK INTEGRATOR CONTROL IS ESSENTIAL. IN C THIS CASE, BRKF45 ATTEMPTS TO ADVANCE THE SOLUTION A C SINGLE BLOCK IN THE DIRECTION OF TEND EACH TIME IT IS C CALLED. SINCE THIS MODE OF OPERATION RESULTS IN EXTRA C COMPUTING OVERHEAD, IT SHOULD BE AVOIDED UNLESS NEEDED. C C C*********************************************************************** C OUTPUT FROM BRKF45 C*********************************************************************** C C Y(*) -- COMPUTED SOLUTION APPROXIMATION AT T. C T -- VALUE OF THE INDEPENDENT VARIABLE WHERE THE SOLUTION IS C REPORTED. C IFLAG = 2 -- SUCCESSFUL RETURN. EITHER THE INTEGRATION REACHED C T=TEND OR A SUCCESSFUL INTERPOLATION HAS BEEN C PERFORMED AT T=TOUT. IF T.EQ.TEND THEN THE C INTEGRATION IS FINISHED. IF NOT, THE CODE SHOULD BE C CALLED WITH THE NEXT VALUE OF TOUT AND WITH IFLAG=+2 C FOR NORMAL INTEGRATION OR IFLAG=-2 FOR ONE-BLOCK C INTEGRATION. C =-2 -- A SINGLE SUCCESSFUL BLOCK IN THE DIRECTION OF TEND C HAS BEEN TAKEN. NORMAL MODE FOR CONTINUING C INTEGRATION ONE BLOCK AT A TIME. C = 3 -- INTEGRATION WAS NOT COMPLETED BECAUSE RELATIVE ERROR C TOLERANCE WAS TOO SMALL. RELERR HAS BEEN INCREASED C APPROPRIATELY FOR CONTINUING. C = 4 -- INTEGRATION WAS NOT COMPLETED BECAUSE MORE THAN C 18000 DERIVATIVE EVALUATIONS WERE NEEDED. THIS C IS APPROXIMATELY 2000 BLOCKS. C = 5 -- INTEGRATION WAS NOT COMPLETED BECAUSE SOLUTION C VANISHED MAKING A PURE RELATIVE ERROR TEST C IMPOSSIBLE. MUST USE NON-ZERO ABSERR TO CONTINUE. C USING THE ONE-BLOCK INTEGRATION MODE FOR ONE BLOCK C IS A GOOD WAY TO PROCEED. C = 6 -- INTEGRATION WAS NOT COMPLETED BECAUSE REQUESTED C ACCURACY COULD NOT BE ACHIEVED USING SMALLEST C ALLOWABLE STEPSIZE. USER MUST INCREASE THE ERROR C TOLERANCE BEFORE CONTINUED INTEGRATION CAN BE C ATTEMPTED. C = 7 -- INVALID INPUT PARAMETERS C THIS INDICATOR OCCURS IF ANY OF THE FOLLOWING IS C SATISFIED - NEQN .LE. 0 C T=TEND C RELERR OR ABSERR .LT. 0. C ABS(IFLAG) .LT. 1 OR .GT.7 C WORK(*),IWORK(*) -- INFORMATION WHICH IS USUALLY OF NO INTEREST C TO THE USER BUT NECESSARY FOR SUBSEQUENT CALLS. C WORK(1),...,WORK(NEQN) CONTAIN THE SOLUTION VECTOR C AND WORK(NEQN+1),...,WORK(2*NEQN) CONTAIN THE C DERIVATIVE VECTOR AT THE END POINT (WHICH IS ITSELF C CONTAINED IN WORK(2*NEQN+1)) OF THE BLOCK STEP JUST C COMPUTED. WORK(2*NEQN+2) CONTAINS THE STEPSIZE H C JUST USED. (THIS IS THE STEPSIZE BEING USED BY THE C INTERPOLANT OVER THIS BLOCK.) WORK(2*NEQN+3) C CONTAINS THE STEPSIZE H TO BE ATTEMPTED ON THE NEXT C BLOCK. IWORK(1) CONTAINS THE DERIVATIVE EVALUATION C COUNTER. C C C*********************************************************************** C SUBSEQUENT CALLS TO BRKF45 C*********************************************************************** C C SUBROUTINE BRKF45 RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE C THE INTEGRATION. AFTER THE CODE REPORTS A SUCCESSFUL SOLUTION AT C TOUT (INDICATED BY IFLAG=2), THE USER NEEDS TO DEFINE A NEW TOUT C BEFORE SIMPLY CALLING BRKF45 AGAIN TO CONTINUE IN THE NORMAL MODE. C (BUT THE USER MUST FIRST RESET IFLAG TO -2 TO CONTINUE IN THE C ONE-BLOCK INTEGRATOR MODE.) C IF THE INTEGRATION WAS NOT COMPLETED BUT THE USER STILL WANTS TO C CONTINUE (IFLAG=3,4), HE JUST CALLS BRKF45 AGAIN. IN THE CASE C IFLAG=3 THE RELERR PARAMETER HAS BEEN ADJUSTED APPROPRIATELY FOR C CONTINUING THE INTEGRATION. IN THE CASE OF IFLAG=4 THE FUNCTION C COUNTER WILL BE RESET TO 0 AND ANOTHER 18000 FUNCTION EVALUATIONS C ARE ALLOWED. C HOWEVER,IN THE CASE IFLAG=5, THE USER MUST FIRST ALTER THE ERROR C CRITERION TO USE A POSITIVE VALUE OF ABSERR BEFORE INTEGRATION CAN C PROCEED. IF HE DOES NOT,EXECUTION IS TERMINATED. C ALSO,IN THE CASE IFLAG=6, IT IS NECESSARY FOR THE USER TO RESET C IFLAG TO 2 (OR -2 WHEN THE ONE-BLOCK INTEGRATION MODE IS BEING C USED) AS WELL AS INCREASING EITHER ABSERR,RELERR OR BOTH BEFORE C THE INTEGRATION CAN BE CONTINUED. IF THIS IS NOT DONE, EXECUTION C WILL BE TERMINATED. THE OCCURRENCE OF IFLAG=6 INDICATES A TROUBLE C SPOT(SOLUTION IS CHANGING RAPIDLY,SINGULARITY MAY BE PRESENT) AND C IT OFTEN IS INADVISABLE TO CONTINUE. C IF IFLAG=7 IS OBTAINED, INTEGRATION CAN NOT BE CONTINUED UNLESS C THE INVALID INPUT PARAMETERS ARE CORRECTED. C IT SHOULD BE NOTED THAT THE ARRAYS WORK,IWORK CONTAIN INFORMATION C REQUIRED FOR SUBSEQUENT INTEGRATION. ACCORDINGLY, WORK AND IWORK C SHOULD NOT BE ALTERED. C C C*********************************************************************** C USER CALLS TO THE INTERPOLANT ROUTINE EXTRA C*********************************************************************** C C SUBROUTINE EXTRA CAN ALSO BE CALLED BY THE USER, IN CONJUNCTION C WITH USAGE OF BRKF45, TO PROVIDE APPROXIMATE SOLUTIONS AT 'OFF-STEP C POINTS' BY USE OF THE INTERPOLATING POLYNOMIAL. WHILE BRKF45 C HANDLES THE USUAL SITUATIONS, IT CAN BE HELPFUL TO THE USER TO BE C ABLE TO ACCESS THE INTERPOLANT DIRECTLY, SUCH AS WHEN DOING ROOT C FINDING. ALSO, IT IS POSSIBLE THAT THE USER MAY HAVE SOME NEED FOR C EXTRAPOLATING OUTSIDE OF THE BLOCK STEP ON WHICH THE UNDERLYING C INTERPOLANT IS BASED. BRKF45 WILL NOT DO THIS. C THE FORM OF THE USAGE CALL IS C C CALL EXTRA ( NEQN, WORK(1), WORK(NEQN+1), WORK(2*NEQN+1), C WORK(2*NEQN+2), WORK(2*NEQN+4), WORK(3*NEQN+4), C WORK(4*NEQN+4), WORK(5*NEQN+4), TEX, YEX, YPEX ) C C WHERE YEX AND YPEX ARE THE SOLUTION VECTOR AND DERIVATIVE VECTOR C APPROXIMATIONS DEFINED BY THE INTERPOLANT AT THE POINT TEX, AND C WORK IS THE WORKING ARRAY SET UP BY BRKF45. C C C*********************************************************************** C MORE ABOUT THE INTERPOLANT ROUTINE EXTRA C*********************************************************************** C C THE FIRST DERIVATIVE APPROXIMATIONS GIVEN BY THE INTERPOLANT C CAN BE MADE LESS SENSITIVE TO ROUNDING ERRORS. IF VERY ACCURATE C FIRST DERIVATIVE APPROXIMATIONS ARE REQUIRED, THEN THREE C CHANGES SHOULD BE MADE TO THE CODE. THESE CHANGES ARE CLEARLY C DOCUMENTED IN THE CODE; JUST SEARCH FOR THE WORD 'CHANGES'. C (ONE LINE OF THE SUBROUTINE RKFC MUST BE REPLACED, ONE LINE C MUST BE ADDED TO THE SUBROUTINE RKFC, AND THE SUBROUTINE C EXTRA MUST BE REPLACED.) C FOR AN EXPLANATION OF THESE MODIFICATIONS SEE C 'REMARK ON ALGORITHM 669', BY D. HIGHAM; TO APPEAR IN TOMS. C C IF THE MODIFICATIONS ARE MADE, THEN THE CALL TO EXTRA DESCRIBED C ABOVE SHOULD BE CHANGED TO C C CALL EXTRA ( NEQN, WORK(6*NEQN+4), WORK(NEQN+1), WORK(2*NEQN+1), C WORK(2*NEQN+2), WORK(2*NEQN+4), WORK(3*NEQN+4), C WORK(4*NEQN+4), WORK(5*NEQN+4), TEX, YEX, YPEX ) C C C*********************************************************************** C LOGICAL ENDPNT,BLKOUT DIMENSION Y(NEQN),YP(NEQN),WORK(*),IWORK(5) EXTERNAL FCN C C COMPUTE INDICES FOR THE SPLITTING OF THE WORK ARRAY C KW = 1 KWP = KW + NEQN KX = KWP + NEQN KHI = KX + 1 KH = KHI + 1 KY1 = KH + 1 KY2 = KY1 + NEQN KF1 = KY2 + NEQN KF2 = KF1 + NEQN KF3 = KF2 + NEQN KF4 = KF3 + NEQN KF7 = KF4 + NEQN KSR = KF7 + NEQN KSA = KSR + 1 KT = KSA + 1 C K ? = KT + 1 C THE WORK SPACE TOTALS 6 + 9*NEQN . C IF (ABS(IFLAG) .NE. 1) THEN ENDPNT = (IWORK(4) .EQ. -1) BLKOUT = (IWORK(5) .EQ. -1) END IF C C*********************************************************************** C THIS INTERFACING ROUTINE MERELY RELIEVES THE USER OF A LONG C CALLING LIST VIA THE SPLITTING APART OF TWO WORKING STORAGE C ARRAYS. IF THIS IS NOT COMPATIBLE WITH THE USERS COMPILER, C HE MUST USE RKFC DIRECTLY. C*********************************************************************** C CALL RKFC (FCN,NEQN,Y,T,TOUT,TEND,YP,RELERR,ABSERR,IFLAG,WORK(KW), * WORK(KWP),WORK(KX),WORK(KHI),WORK(KH),WORK(KY1),WORK(KY2), * WORK(KF1),WORK(KF2),WORK(KF3),WORK(KF4),WORK(KF7),WORK(KSR), * WORK(KSA),WORK(KT),IWORK(1),IWORK(2),IWORK(3),ENDPNT,BLKOUT) C IWORK(4) = 0 IF (ENDPNT) IWORK(4) = -1 IWORK(5) = 0 IF (BLKOUT) IWORK(5) = -1 C RETURN END SUBROUTINE RKFC (FCN,NEQN,Y,T,TOUT,TEND,YP,RELERR,ABSERR,IFLAG, * W,WP,X,HINT,H,Y1,Y2,F1,F2,F3,F4,F7,SAVRE,SAVAE, * TOLD,NFE,JFLAG,KFLAG,ENDPNT,BLKOUT) C C TWO STEP BLOCK RUNGE-KUTTA FEHLBERG METHOD. C A STANDARD 5(4) FORMULA IS USED AT THE FIRST POINT IN THE BLOCK C AND A 6(4) FORMULA IS USED AT THE SECOND POINT. C LOGICAL HFAILD,ENDPNT,INTERP,BLKOUT DIMENSION Y(NEQN),YP(NEQN),Y1(NEQN),Y2(NEQN),F1(NEQN),F2(NEQN), * F3(NEQN),F4(NEQN),F7(NEQN),W(NEQN),WP(NEQN) EXTERNAL FCN C C*********************************************************************** C C COEFFICIENTS DEFINING THE METHOD ARE SET IN THE FOLLOWING C PARAMETER STATEMENTS. C PARAMETER( C2=0.25E+0, C3=3.0E+0/8.0E+0, C4=12.0E+0/13.0E+0, * C6=0.5E+0, C8=3.0E+0/2.0E+0, C9=2.0E+0 ) PARAMETER( A21=1.0E+0/4.0E+0) PARAMETER (A31=3.0E+0/32.0E+0, A32=9.0E+0/32.0E+0 ) PARAMETER( A41=1932.0E+0/2197.0E+0, A42=-7200.0E+0/2197.0E+0, * A43=7296.0E+0/2197.0E+0 ) PARAMETER( A51=439.0E+0/216.0E+0, A52=-8.0E+0, * A53=3680.0E+0/513.0E+0, A54=-845.0E+0/4104.0E+0) PARAMETER( A61=-8.0E+0/27.0E+0, A62=2.0E+0, * A63=-3544.0E+0/2565.0E+0, A64=1859.0E+0/4104.0E+0, * A65=-11.0E+0/40.0E+0) PARAMETER( B1B=16.0E+0/135.0E+0, B3B=6656.0E+0/12825.0E+0, * B4B=28561.0E+0/56430.0E+0, B5B=-9.0E+0/50.0E+0, * B6B=2.0E+0/55.0E+0) PARAMETER( ERC11=1.0E+0/360.0E+0, ERC13=-128.0E+0/4275.0E+0, * ERC14=-2197.0E+0/75240.0E+0, ERC15=1.0E+0/50.0E+0) C C THE ABOVE DEFINE THE COEFFICIENTS FOR BRKF45 USED TO GENERATE C THE SOLUTION AT THE FIRST BLOCK POINT. BELOW ARE ADDITIONAL C COEFFICIENTS NEEDED TO GENERATE THE SOLUTION AT THE SECOND BLOCK C POINT. C PARAMETER( B1=931.0E+0/6480.0E+0, B3=315392.0E+0/1500525.0E+0, * B4=371293.0E+0/615600.0E+0, B5=1.0E+0/50.0E+0, * B6=0.4E+0, B7=-4.0E+0/15.0E+0, * B8=85006.0E+0/115425.0E+0, B9=239.0E+0/1560.0E+0) PARAMETER( A81=-119397029895.0E+0/151948225000.0E+0, * A82=78390.0E+0/29081.0E+0, * A83=-51517464.0E+0/132821875.0E+0, * A84=-3780749193.0E+0/1168832500.0E+0, * A85=79268193.0E+0/55925000.0E+0, * A86=-11370591.0E+0/15379375.0E+0, * A87=5670.0E+0/2237.0E+0) PARAMETER( A91=23406188597.0E+0/8429231250.0E+0, * A92=-62928.0E+0/13623.0E+0, * A93=-31066887488.0E+0/5747203125.0E+0, * A94=164486461399.0E+0/8429231250.0E+0, * A95=-70336084.0E+0/11203125.0E+0, * A96=185680664.0E+0/24646875.0E+0, * A97=-3385330161.0E+0/243117160.0E+0, * A98=232648.0E+0/96795.0E+0) C C NEXT WE DEFINE COEFFICIENTS FOR THE ERROR ESTIMATE FORMULA. C PARAMETER(ERC21=0.22490537549467E-1,ERC23=-0.18660147916167E+0, * ERC24=0.224620349650850E+0, ERC25=-0.38462277720213E-1, * ERC26=9.0E+0/50.0E+0, ERC27=-7.0E+0/30.0E+0, * ERC28=0.36286203014893E-1, ERC29=-0.005E+0) C C IN WHAT FOLLOWS WE GIVE THE RATIONAL FORMS OF THESE DECIMAL C CONSTANTS THAT DEFINE THE ERROR COEFFICIENTS. C C B3E = 1067091077380.0E+0/1829119027671.0E+0 C B4E = 3284168845918.0E+0/21339721989495.0E+0 C B5E = 110317750789.0E+0/240996319200.0E+0- C * 4448925830089.0E+0/12329617149531.0E+0 C B6E = 1.0E+0/25.0E+0 C B7E = 0.2E+0 C B8E = 239992027043.0E+0/361494478800.0E+0 C B9E = 1273.0E+0/7800.0E+0 C B1E = 2.0E+0-B3E-B4E-B5E-B6E-B7E-B8E-B9E C ERC21 = (B1-B1E)/2.0E+0 C ERC23 = (B3-B3E)/2.0E+0 C ERC24 = (B4-B4E)/2.0E+0 C ERC25 = (B5-B5E)/2.0E+0 C ERC26 = 9.0E+0/50.0E+0 C ERC27 = -7.0E+0/30.0E+0 C ERC28 = (B8 - B8E)/2.0E+0 C ERC29 = -0.005E+0 C C IN WHAT FOLLOWS WE GIVE DECIMAL APPROXIMATIONS TO THE C COEFFICIENTS NEEDED TO COMPUTE THE SOLUTION OVER THE SECOND C STEP OF THE BLOCK. C C B1=0.14367283950617 C B3=0.21018776761467 C B4=0.60314002599090 C B5=0.02 C B6=0.4 C B7=-0.26666666666667 C B8=0.73646090534979 C B9=0.15320512820513 C A81=-0.78577443004023 C A82=2.6955744300402 C A83=-0.38786882055384 C A84=-3.2346372923408 C A85=1.4174017523469 C A86=-0.73934025277360 C A87=2.5346446133214 C A91=2.7767880489695 C A92=-4.6192468619247 C A93=-5.4055662923868 C A94=19.513815260318 C A95=-6.2782557545328 C A96=7.5336392138963 C A97=-13.924686192451 C A98=2.4035125781290 C B1E=0.098691764407238 C B3E=0.58339072593800 C B4E=0.1538993266892 C B5E=0.096924555440426 C B6E=0.04 C B7E=0.2 C B8E=0.66388849932001 C B9E=0.16320512820513 C C THE COEFFICIENTS OF THE RUNGE-KUTTA FORMULA ARE NOW SET. C C*********************************************************************** C C THE COMPUTER UNIT ROUNDOFF ERROR U IS THE SMALLEST POSITIVE VALUE C REPRESENTABLE IN THE MACHINE SUCH THAT 1.+ U .GT. 1. C VALUES TO BE USED ARE C U = 9.5E-7 FOR IBM 360/370 C U = 1.2E-7 FOR DEC VAX (SINGLE PRECISION) C U = 1.5E-8 FOR UNIVAC 1108 C U = 7.5E-9 FOR PDP-10 C U = 7.1E-15 FOR CDC 6000 SERIES C U = 2.2D-16 FOR IBM 360/370 (DOUBLE PRECISION) C U = 2.8D-17 FOR DEC VAX (DOUBLE PRECISION) C C DATA U / 1.2E-7 / DATA U / 7.1E-15 / C C*********************************************************************** C C REMIN IS A TOLERANCE THRESHOLD WHICH IS ALSO DETERMINED BY THE C INTEGRATION METHOD. IN PARTICULAR, A FIFTH ORDER METHOD WILL C GENERALLY NOT BE CAPABLE OF DELIVERING ACCURACIES NEAR LIMITING C PRECISION ON COMPUTERS WITH LONG WORDLENGTHS. THIS DOES NOT HAVE TO C BE CHANGED FOR DIFFERENT MACHINES. C DATA REMIN / 1.E-12 / C C*********************************************************************** C C THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER C OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE. C DATA MAXNFE / 18000 / C C*********************************************************************** C C C CHECK INPUT PARAMETERS C C IF (NEQN.LT.1) GO TO 10 IF ((RELERR.LT.0.0E+0).OR.(ABSERR.LT.0.0E+0)) GO TO 10 IF (T.EQ.TEND) GO TO 10 MFLAG = ABS(IFLAG) IF (MFLAG .EQ. 1) THEN TOLD = T IF (TOUT .EQ. T) GO TO 20 IF (SIGN(1.0E+0,TEND-T) .NE. SIGN(1.0E+0,TOUT-T)) GO TO 10 ELSE IF (TOUT .EQ. TOLD) GO TO 20 IF (SIGN(1.0E+0,TEND-TOLD) .NE. SIGN(1.0E+0,TOUT-TOLD)) * GO TO 10 END IF IF ((MFLAG.GE.1).AND.(MFLAG.LE.7)) GO TO 20 C C INVALID INPUT, RETURN C 10 IFLAG = 7 RETURN C C IS THIS THE FIRST CALL? IF SO JUMP TO 70. C 20 IF (MFLAG.EQ.1) GO TO 70 C C CHECK CONTINUATION POSSIBILITIES C IF (MFLAG.NE.2) GO TO 30 C C IFLAG = +2 OR -2 C IF (KFLAG.LT.3) GO TO 70 IF (KFLAG.EQ.3) GO TO 60 IF (KFLAG.EQ.4) GO TO 50 IF ((KFLAG.EQ.5).AND.(ABSERR.EQ.0.0E+0)) GO TO 40 IF ((KFLAG.EQ.6).AND.(RELERR.LE.SAVRE).AND.(ABSERR.LE.SAVAE)) * GO TO 40 GO TO 70 C C IFLAG = 3,4,5,6 OR 7 C 30 IF (IFLAG.EQ.3) GO TO 60 IF (IFLAG.EQ.4) GO TO 50 IF ((IFLAG.EQ.5).AND.(ABSERR.GT.0.0E+0)) GO TO 60 C C INTEGRATION CANNOT BE CONTINUED SINCE USER DID NOT RESPOND TO C THE INSTRUCTIONS PERTAINING TO IFLAG=5,6 OR 7 C 40 STOP C C C*********************************************************************** C C RESET FUNCTION EVALUATION COUNTER C 50 NFE = 0 IF (MFLAG.EQ.2) GO TO 70 C C RESET FLAG VALUE FROM PREVIOUS CALL C 60 IFLAG = JFLAG IF (KFLAG.EQ.3) MFLAG = ABS(IFLAG) C C SAVE INPUT IFLAG AND SET CONTINUATION FLAG VALUE FOR SUBSEQUENT C INPUT CHECKING C 70 JFLAG = IFLAG KFLAG = 0 C C SAVE RELERR AND ABSERR FOR CHECKING INPUT ON SUBSEQUENT CALLS C SAVRE = RELERR SAVAE = ABSERR C C RESTRICT RELATIVE ERROR TOLERANCE TO BE AT LEAST AS LARGE AS C 2U+REMIN TO AVOID LIMITING PRECISION DIFFICULTIES ARISING FROM C IMPOSSIBLE ACCURACY REQUESTS. IF TOLERANCE TOO SMALL, INCREASE C AND RETURN. C RER = 2.0E+0*U+REMIN IF (RELERR .LT. RER) THEN C C RELATIVE ERROR TOLERANCE TOO SMALL C RELERR = RER IFLAG = 3 KFLAG = 3 RETURN END IF C U26 = 26.0E+0*U IF (MFLAG.NE.1) GO TO 100 C C*********************************************************************** C C INITIALIZATION -- C DEFINE INTEGRATION INDEPENDENT VARIABLE X C EVALUATE INITIAL DERIVATIVES C SET UP WORKING ARRAYS FOR INTEGRATION VARIABLES C SET COUNTER FOR FUNCTION EVALUATIONS,NFE C ESTIMATE STARTING STEPSIZE C X = T ENDPNT = .FALSE. BLKOUT = .FALSE. C A = T CALL FCN (A,Y,YP) NFE = 1 DO 80 N = 1, NEQN W(N) = Y(N) WP(N) = YP(N) 80 CONTINUE C C COMPUTE INITIAL STEPLENGTH. C DT = TOUT - T IF (DT .EQ. 0.0E+0) DT = TEND - T H = ABS(DT) TOLN = 0.0E+0 DO 90 K = 1, NEQN TOL = RELERR*ABS(Y(K))+ABSERR IF (TOL.LE.0.0E+0) GO TO 90 TOLN = TOL YPK = ABS(YP(K)) IF (YPK*H**5.GT.TOL) H = (TOL/YPK)**0.2E+0 90 CONTINUE IF (TOLN.LE.0.0E+0) H = 0.0E+0 H = MAX(H,U26*MAX(ABS(T),ABS(DT))) JFLAG = SIGN(2,IFLAG) H = SIGN(H,DT) C C INITIAL STEPLENGTH NOW COMPUTED. COMPUTE FIRST SOLUTION. C C*********************************************************************** C 100 CONTINUE C C TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION, C SCALE THE ERROR TOLERANCES. C SCALE = 2.0E+0/RELERR AE = SCALE*ABSERR C C SET SAFETY FACTOR FOR STEPSIZE ADJUSTMENT, BASED ON TOLERANCES. C TOLER = MAX(ABSERR,RELERR) SF = 0.85E+0 IF (TOLER .GE. 1.E-5) SF = 0.8E+0 IF (TOLER .LE. 1.E-9) SF = 0.9E+0 C C RESTORE INTEGRATION VARIABLE TO END OF LAST BLOCK STEP TAKEN C AND SET THE SIGN OF THE DIRECTION OF INTEGRATION. C T = X DTSIGN = SIGN(1.0E+0,TEND-T) C C HAVE WE ALREADY INTEGRATED PAST THE PRESENT DATA OUTPUT POINT? C IF SO JUMP TO 390 AND PERFORM INTERPOLATION. C IF NOT, SEE IF WE HAVE REACHED THE END POINT OF INTEGRATION. C IF NOT, SEE IF RESULTS AT THE END OF THE BLOCK STEP NEED TO BE C REPORTED. C IF ((TOUT-T)*DTSIGN .LE. 0.0E+0) THEN IF (TOUT .EQ. T) THEN IFLAG = 2 ENDPNT = .FALSE. GO TO 400 ELSE GO TO 390 END IF END IF IF (ENDPNT) THEN IFLAG = 2 ENDPNT = .FALSE. GO TO 400 END IF IF (IFLAG .EQ. -2 .AND. BLKOUT) THEN BLKOUT = .FALSE. GO TO 400 END IF C C*********************************************************************** C*********************************************************************** C BLOCK BY BLOCK INTEGRATION C 150 CONTINUE C C SEE IF WE ARE TOO CLOSE TO THE END POINT. IF SO, DO LINEAR C EXTRAPOLATION AND RETURN. C IF (ABS(TEND-T) .LE. U26*ABS(T)) THEN IF (DTSIGN*(TEND-TOUT) .GT. 0.0E+0) THEN DT = TOUT - T ENDPNT = .FALSE. T = TOUT ELSE DT = TEND - T ENDPNT = .TRUE. T = TEND X = TEND END IF DO 160 K = 1, NEQN Y(K) = W(K)+DT*WP(K) 160 CONTINUE CALL FCN (T,Y,YP) NFE = NFE+1 IF (ENDPNT) THEN DO 170 N=1,NEQN W(N) = Y(N) WP(N) = YP(N) 170 CONTINUE ENDPNT = .FALSE. END IF IFLAG = 2 RETURN END IF C C C SET SMALLEST ALLOWABLE STEPSIZE AND STEP FAILURE FLAG. C ADJUST STEPSIZE IF NECESSARY TO HIT THE END POINT OF INTEGRATION. C HMIN = U26*ABS(T) HFAILD = .FALSE. HSTOP = 0.5E+0*(TEND-T) IF (ABS(HSTOP) .LE. ABS(H)) THEN ENDPNT = .TRUE. H = HSTOP END IF C C*********************************************************************** C CORE INTEGRATOR FOR A SINGLE BLOCK C*********************************************************************** C THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW IN C COMPUTING THE ERROR TOLERANCE FUNCTION ERRTOL. C TO AVOID PROBLEMS WITH ZERO CROSSINGS,RELATIVE ERROR IS MEASURED C USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION AT THE C BEGINNING AND END POINTS OF A BLOCK. C TO DISTINGUISH THE VARIOUS ARGUMENTS, H IS NOT PERMITTED C TO BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T. C PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED TO C SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID EXCESSIVE C CHATTERING ON PROBLEMS HAVING DISCONTINUITIES. C*********************************************************************** C C C TEST NUMBER OF DERIVATIVE FUNCTION EVALUATIONS. C IF OKAY,TRY TO ADVANCE THE INTEGRATION FROM T TO T+H C 200 IF (NFE .GT. MAXNFE) THEN C C TOO MUCH WORK C IFLAG = 4 KFLAG = 4 GO TO 400 END IF C C ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H C DO 210 N = 1, NEQN Y1(N) = W(N)+A21*H*WP(N) 210 CONTINUE CALL FCN (T+C2*H,Y1,F2) DO 220 N = 1, NEQN Y1(N) = W(N)+H*(A31*WP(N)+A32*F2(N)) 220 CONTINUE CALL FCN (T+C3*H,Y1,F3) DO 230 N = 1, NEQN Y1(N) = W(N)+H*(A41*WP(N)+A42*F2(N)+A43*F3(N)) 230 CONTINUE CALL FCN (T+C4*H,Y1,F4) DO 240 N = 1, NEQN Y1(N) = W(N)+H*(A51*WP(N)+A52*F2(N)+A53*F3(N)+A54*F4(N)) 240 CONTINUE CALL FCN (T+H,Y1,Y) DO 250 N = 1, NEQN Y1(N) = W(N)+H*(A61*WP(N)+A62*F2(N)+A63*F3(N)+A64*F4(N)+ * A65*Y(N)) 250 CONTINUE CALL FCN (T+C6*H,Y1,YP) NFE = NFE+5 EEOET = 0.0E+0 DO 270 N = 1, NEQN Y1(N) = W(N)+H*(B1B*WP(N)+B3B*F3(N)+B4B*F4(N)+B5B*Y(N)+ * B6B*YP(N)) ERRTOL = ABS(W(N))+ABS(Y1(N))+AE IF (ERRTOL .LE. 0.0E+0) THEN C C INAPPROPRIATE ERROR TOLERANCE C IFLAG = 5 KFLAG = 5 GO TO 400 END IF EZ = ABS(H*(ERC11*WP(N)+ERC13*F3(N)+ERC14*F4(N)+ERC15*Y(N)+ * B6B*YP(N))) EEOET = MAX(EEOET,EZ/ERRTOL) 270 CONTINUE ESTTOL = EEOET*SCALE C C CHECK THE ERROR ESTIMATE. IF STEPLENGTH HAS FAILED FOR FIRST STEP C IN THE BLOCK, GO AND COMPUTE A SMALLER STEP FOR A RE-TRY. C IF (ESTTOL .GT. 1.0E+0) GO TO 320 C C C INTEGRATION SUCCESSFUL FOR FIRST STEP. NOW INTEGRATE OVER C THE SECOND STEP IN THE BLOCK C CALL FCN (T+H,Y1,F1) DO 290 N = 1, NEQN Y2(N) = W(N)+H*(A81*WP(N)+A82*F2(N)+A83*F3(N)+A84*F4(N)+ * A85*Y(N)+A86*YP(N)+A87*F1(N)) 290 CONTINUE CALL FCN (T+C8*H,Y2,F7) DO 300 N = 1, NEQN Y2(N) = W(N)+H*(A91*WP(N)+A92*F2(N)+A93*F3(N)+A94*F4(N)+ * A95*Y(N)+A96*YP(N)+A97*F1(N)+A98*F7(N)) 300 CONTINUE CALL FCN (T+C9*H,Y2,F2) NFE = NFE+3 EEOET = 0.0E+0 DO 310 N = 1, NEQN EE = ABS(H*(ERC21*WP(N)+ERC23*F3(N)+ERC24*F4(N)+ERC25*Y(N)+ * ERC26*YP(N)+ERC27*F1(N)+ERC28*F7(N)+ERC29*F2(N))) Y2(N) = W(N)+H*(B1*WP(N)+B3*F3(N)+B4*F4(N)+B5*Y(N)+B6*YP(N)+ * B7*F1(N)+B8*F7(N)+B9*F2(N)) ERRTOL = ABS(Y2(N))+ABS(Y1(N))+AE IF (ERRTOL .LE. 0.0E+0) THEN C C INAPPROPRIATE ERROR TOLERANCE C IFLAG = 5 KFLAG = 5 GO TO 400 END IF EEOET = MAX(EEOET,EE/ERRTOL) 310 CONTINUE ESTTOL = EEOET*SCALE C C THE FIRST OF THE THREE CHANGES NEEDED TO MAKE THE INTERPOLANT C DERIVATIVE MORE STABLE IS TO UNCOMMENT THE FOLLOWING DO-LOOP C C DO 315 N = 1, NEQN C Y1(N) = B1B*WP(N)+B3B*F3(N)+B4B*F4(N)+B5B*Y(N)+ C * B6B*YP(N) C F3(N) = B1*WP(N)+B3*F3(N)+B4*F4(N)+B5*Y(N)+B6*YP(N)+ C * B7*F1(N)+B8*F7(N)+B9*F2(N) C 315 CONTINUE C C CHECK THE ERROR ESTIMATE OVER THE SECOND STEP OF THE BLOCK. C IF (ESTTOL .LE. 1.0E+0) GO TO 330 C C UNSUCCESSFUL BLOCK C REDUCE THE STEPSIZE , TRY AGAIN C THE DECREASE IS LIMITED TO A FACTOR OF ABOUT 1/10 C 320 HFAILD = .TRUE. S = 0.1E+0 ENDPNT = .FALSE. IF (ESTTOL.LT.1.0E+5) S = SF/ESTTOL**0.2E+0 H = S*H IF (ABS(H) .LT. HMIN) THEN C C REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE C IFLAG = 6 KFLAG = 6 GO TO 400 END IF GO TO 200 C C C SUCCESSFUL BLOCK. CHECK FOR NEED TO INTERPOLATE. C STORE SOLUTION AT T+2*H. C 330 T = T+2.0E+0*H IF (ENDPNT) T = TEND TOLD = X X = T HINT = H INTERP = .FALSE. IF ((T-TOUT)*DTSIGN .GE. 0.0E+0) INTERP = .TRUE. IF (INTERP) THEN DO 340 N=1,NEQN SWAP = W(N) W(N) = Y2(N) Y2(N) = SWAP F2(N) = WP(N) 340 CONTINUE ELSE DO 350 N=1,NEQN W(N) = Y2(N) 350 CONTINUE END IF C A = T CALL FCN (A,W,WP) NFE = NFE+1 C C CHOOSE NEXT STEPSIZE C THE INCREASE IS LIMITED TO A FACTOR OF ABOUT 10 C IF STEP FAILURE HAS JUST OCCURED, NEXT STEP IS NOT ALLOWED C TO INCREASE. C S = 10.0E+0 IF (ESTTOL.GT.1.E-5) S = SF/ESTTOL**0.2E+0 IF (HFAILD) S = MIN(S,1.0E+0) H = SIGN(MAX(S*ABS(H),HMIN),H) C C HAVE WE INTEGRATED PAST AN OUTPUT POINT? C IF SO, CALL THE INTERPOLATION ROUTINE AT T=TOUT. C IF NOT, SEE IF WE ARE AT THE END POINT OF INTEGRATION. C OTHERWISE, CHECK IF USER WANTS SOLUTIONS AT THE END OF THE BLOCK C STEP, OR ELSE CONTINUE THE INTEGRATION. C IF (INTERP) GO TO 390 IF (ENDPNT) THEN IFLAG = 2 ENDPNT = .FALSE. GO TO 400 END IF IF (IFLAG .LT. 0) THEN IFLAG = -2 GO TO 400 ELSE GO TO 150 END IF C C INTERPOLATE TO GET DATA AT OFF-STEP POINT AND RETURN. C C THE SECOND OF THE THREE CHANGES NEEDED TO MAKE THE INTERPOLANT C DERIVATIVE MORE STABLE IS TO REPLACE STATEMENT 390 BY THE C VERSION IN THE COMMENT STATEMENT BELOW C 390 CALL EXTRA (NEQN,W,WP,X,HINT,Y1,Y2,F1,F2,TOUT,Y,YP) C 390 CALL EXTRA (NEQN,F3,WP,X,HINT,Y1,Y2,F1,F2,TOUT,Y,YP) T = TOUT IF (IFLAG .LT. 0 .AND. TOUT .NE. X) BLKOUT = .TRUE. IFLAG = 2 RETURN C C RETURN WITH THE SOLUTION AT THE END OF THE LAST SUCCESSFUL C BLOCK STEP. C 400 DO 410 N=1,NEQN Y(N) = W(N) YP(N) = WP(N) 410 CONTINUE RETURN C END C C THE THIRD OF THE THREE CHANGES NEEDED TO MAKE THE INTERPOLANT C DERIVATIVE MORE STABLE IS TO REPLACE THE SUBROUTINE EXTRA BY C THE VERSION IN COMMENT STATEMENTS BELOW C SUBROUTINE EXTRA (NEQN,Y2,F2,T,H,Y1,Y,F1,F,TEX,YEX,YPEX) C C T IS THE INDEPENDENT VARIABLE; H IS THE STEP SPACING FOR VALUES C Y,Y1,Y2; TEX IS THE POINT WHERE THE OUTPUT IS REQUIRED; C F,F1 AND F2 ARE THE DERIVATIVES AT T-2H, T-H AND T. C YEX WILL HOLD THE APPROXIMATE SOLUTION AT TEX AND YPEX WILL C HOLD THE DERIVATIVE APPROXIMATION AT THIS POINT. C DIMENSION Y(NEQN),Y1(NEQN),Y2(NEQN),F(NEQN),F1(NEQN),F2(NEQN), * YEX(NEQN),YPEX(NEQN) C C PERFORM QUINTIC INTERPOLATION BASED ON THE VALUES C Y(N),Y1(N),Y2(N),F(N),F1(N),F2(N) AT THE POINTS C T-2H, T-H, T. THIS POLYNOMIAL IS EVALUATED AT THE POINT TEX C TO OBTAIN THE APPROXIMATE VALUE OF Y(TEX) AND THIS IS STORED C IN THE ARRAY YEX. THE DERIVATIVE OF THE INTERPOLATING POLYNOMIAL C IS ALSO EVALUATED AT TEX TO OBTAIN AN APPROXIMATION TO DY/DT C AT THIS POINT AND THIS APPROXIMATION IS STORED IN YPEX. C C SIG = (TEX-T)/H+2.0E+0 DO 10 N = 1, NEQN HF = H*F(N) HF1 = H*F1(N) HF2 = H*F2(N) A1 = 2.0E+0*(HF2+HF)+6.0E+0*(Y(N)-Y2(N))+8.0E+0*HF1 A2 = Y2(N)+4.0E+0*(Y1(N)-HF1)-5.0E+0*Y(N)-2.0E+0*HF A3 = HF1+HF+2.0E+0*(Y(N)-Y1(N)) A4 = Y1(N)-Y(N)-HF YEX(N) = ((((A1*0.125E+0*(SIG-2.0E+0)+0.25E+0*A2)*(SIG-1.0E+0)+ * A3)*(SIG-1.0E+0)+A4)*SIG+HF)*SIG+Y(N) YPEX(N) = ((((0.625E+0*A1*SIG+(A2-2.0E+0*A1))*SIG+(1.875E+0*A1- * 1.5E+0*A2+3.0E+0*A3))*SIG+(0.5E+0*(A2-A1)- * 2.0E+0*(A3-A4)))*SIG+HF)/H 10 CONTINUE RETURN END C C C C SUBROUTINE EXTRA (NEQN,YINC2,F2,T,H,YINC1,Y,F1,F,TEX,YEX,YPEX) CC CC T IS THE INDEPENDENT VARIABLE; CC Y IS THE VALUE AT T-2H; THE VALUES AT T-H AND T ARE CC Y + H*YINC1 AND Y + H*YINC2 RESPECTIVELY; CC F,F1 AND F2 ARE THE DERIVATIVES AT T-2H, T-H AND T. CC TEX IS THE POINT WHERE THE OUTPUT IS REQUIRED; CC YEX WILL HOLD THE APPROXIMATE SOLUTION AT TEX AND YPEX WILL CC HOLD THE DERIVATIVE APPROXIMATION AT THIS POINT. CC C DIMENSION Y(NEQN),YINC1(NEQN),YINC2(NEQN),F(NEQN),F1(NEQN), C * F2(NEQN),YEX(NEQN),YPEX(NEQN) CC CC PERFORM QUINTIC INTERPOLATION BASED ON THE VALUES CC Y(N),Y(N)+H*YINC1(N),Y(N)+H*YINC2(N),F(N),F1(N),F2(N) AT THE CC POINTS T-2H, T-H, T. THIS POLYNOMIAL IS EVALUATED AT THE POINT CC TEX TO OBTAIN THE APPROXIMATE VALUE OF Y(TEX) AND THIS IS STORED CC IN THE ARRAY YEX. THE DERIVATIVE OF THE INTERPOLATING POLYNOMIAL CC IS ALSO EVALUATED AT TEX TO OBTAIN AN APPROXIMATION TO DY/DT CC AT THIS POINT AND THIS APPROXIMATION IS STORED IN YPEX. CC CC C SIG = (TEX-T)/H+2.0E+0 C DO 10 N = 1, NEQN C HF = H*F(N) C HF1 = H*F1(N) C HF2 = H*F2(N) C A1 = 2.0E+0*(HF2+HF)-6.0E+0*H*YINC2(N)+8.0E+0*HF1 C A2 = H*YINC2(N)+4.0E+0*H*YINC1(N)-4.0E+0*HF1-2.0E+0*HF C A3 = HF1+HF-2.0E+0*H*YINC1(N) C A4 = H*YINC1(N)-HF C YEX(N) = ((((A1*0.125E+0*(SIG-2.0E+0)+0.25E+0*A2)*(SIG-1.0E+0)+ C * A3)*(SIG-1.0E+0)+A4)*SIG+HF)*SIG+Y(N) C YPEX(N) = ((((0.625E+0*A1*SIG+(A2-2.0E+0*A1))*SIG+(1.875E+0*A1- C * 1.5E+0*A2+3.0E+0*A3))*SIG+(0.5E+0*(A2-A1)- C * 2.0E+0*(A3-A4)))*SIG+HF)/H C 10 CONTINUE C RETURN C END * C PROGRAM RUN (INPUT,OUTPUT,RESULT,TAPE5=INPUT,TAPE6=RESULT) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C PROGRAM TO INTEGRATE DY/DT=-Y, Y(0)=1 IN THE RANGE 0.LE.T.LE.20 C WITH OUTPUT REQUIRED AT T=1.0,2.0,.....,19.0,20.0. C C THE DIMENSION STATEMENT HAS BEEN SET UP FOR 10 DIFFERENTIAL C EQUATIONS. C DIMENSION Y(10),YP(10),WORK(96),IWORK(5) EXTERNAL FCN C C SET INPUT PARAMETERS TO SUBROUTINE BRKF45 C NEQN = 1 Y(1) = 1.0D+0 T = 0.0D+0 IFLAG = -1 DEL = 1.0D+0 K = 1 TOUT = DFLOAT(K)*DEL TEND = 20.0D+0 RELERR = 1.D-10 ABSERR = 1.D-6 WRITE(6, 50) C C INPUT PARAMETERS SET. C 10 CALL BRKF45 (FCN,NEQN,Y,T,TOUT,TEND,YP,RELERR,ABSERR,IFLAG,WORK, * IWORK) IF (IFLAG.LE.2) GO TO 30 C C IF IFLAG IS GREATER THAN 2 WE HAVE AN INPUT ERROR. C WRITE (6,20) IFLAG,NEQN,IWORK(1),T,TOUT,RELERR,ABSERR 20 FORMAT (1X,'FAILURE',5X,'IFLAG,NEQN,NUM F EVAL = ',2I5,I10/ * 13X,'T,TOUT = ',2D13.4,/13X,'RELERR,ABSERR = ',1P2D10.1) IF (IFLAG .EQ. 3) THEN WRITE (6,25) 25 FORMAT (1X,'RELATIVE ERROR TOLERANCE INCREASED TO CONTINUE') GO TO 10 ELSE STOP END IF 30 CONTINUE IF (T.LT.TOUT) GO TO 10 C C THE INTEGRATION HAS PASSED THE OUTPUT POINT T=TOUT. C A = DEXP(-T) B = -A E1 = Y(1)-A E2 = YP(1)-B C C E1 AND E2 ARE THE ERRORS IN THE INTERPOLATED SOLUTION C AND DERIVATIVE C WRITE (6,60) T,Y(1),YP(1),E1,E2 C C SET NEW VALUE OF THE OUTPUT POINT. C K = K + 1 TOUT = DFLOAT(K)*DEL IF(DABS(TOUT-TEND).LT.ABSERR) TOUT=TEND IF (T.NE.TEND) GO TO 10 50 FORMAT(5X,1HT, 11X, 6H Y, 14X, * 10H DY/DT, 6X, 10HERROR IN Y, 3X, * 14HERROR IN DY/DT) 60 FORMAT (1X,F8.4,1P2D22.10,2D13.2) END SUBROUTINE FCN (T,Y,YP) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION Y(1),YP(1) YP(1) = -Y(1) RETURN END SUBROUTINE BRKF45 (FCN,NEQN,Y,T,TOUT,TEND,YP,RELERR,ABSERR,IFLAG, * WORK,IWORK) C C BLOCK FEHLBERG FOURTH-FIFTH ORDER RUNGE-KUTTA METHOD ADVANCING C A BLOCK OF TWO EQUAL STEPS. THE FORMULA AT THE FIRST STEP IS C 5(4) WHILE AT THE SECOND STEP IT IS 6(4). C C WRITTEN BY J.R.CASH, C DEPARTMENT OF MATHEMATICS, C IMPERIAL COLLEGE, C SOUTH KENSINGTON, LONDON SW7 2AZ, C ENGLAND. C C THIS IS A HEAVILY REVISED VERSION OF RKF45 OF L.F. SHAMPINE AND C H.A. WATTS. C BRKF45 IS PRIMARILY DESIGNED TO SOLVE NON-STIFF AND MILDLY STIFF C INITIAL VALUE ORDINARY DIFFERENTIAL EQUATIONS WHEN DERIVATIVE C EVALUATIONS ARE INEXPENSIVE. BRKF45 USES INTERPOLATION TO PRODUCE C OUTPUT AT 'OFF-STEP POINTS' EFFICIENTLY. BRKF45 SHOULD GENERALLY C NOT BE USED WHEN THE USER IS DEMANDING HIGH ACCURACY. IN SUCH C CASES A GOOD ADAMS CODE WILL OFTEN BE MORE EFFICIENT. C C*********************************************************************** C ABSTRACT C*********************************************************************** C C SUBROUTINE BRKF45 INTEGRATES A SYSTEM OF NEQN FIRST ORDER C ORDINARY DIFFERENTIAL EQUATIONS OF THE FORM C C DY(I)/DT = FCN(T,Y(1),Y(2),...,Y(NEQN)) C C WHERE THE Y(I) ARE KNOWN AT TIME T. C TYPICALLY THE SUBROUTINE IS USED TO INTEGRATE FROM T TO TEND C (WHILE RETURNING ANSWERS AT SPECIFIED OUTPUT POINTS TOUT), BUT C IT CAN ALSO BE USED AS A ONE-BLOCK INTEGRATOR TO ADVANCE THE C SOLUTION A SINGLE BLOCK STEP IN THE DIRECTION OF TEND. ON RETURN C THE PARAMETERS IN THE CALL LIST ARE SET FOR CONTINUING THE C INTEGRATION. THE USER HAS ONLY TO CALL BRKF45 AGAIN (AND PERHAPS C DEFINE A NEW VALUE FOR TOUT). ACTUALLY, BRKF45 IS AN INTERFACING C ROUTINE WHICH CALLS SUBROUTINE RKFC FOR THE SOLUTION. SUBROUTINE C RKFC COMPUTES AN APPROXIMATE SOLUTION OVER ONE BLOCK OF LENGTH 2H. C BRKF45 IS PARTICULARLY USEFUL WHEN OUTPUT IS REQUIRED AT MANY C 'OFF-STEP POINTS' SINCE THE OUTPUT VALUES CAN BE OBTAINED BY C INTERPOLATION. THIS IS IN CONTRAST TO MANY OTHER RUNGE-KUTTA C PROGRAMS WHICH CHOOSE THE STEP SEQUENCE SO AS TO HIT ALL OUTPUT C POINTS EXACTLY AND SO BECOME INEFFICIENT WHEN OUTPUT IS REQUIRED C AT MANY POINTS WITHIN A STEP. C BRKF45 USES THE (5,4); (6,4) BLOCK FORMULA DESCRIBED IN C J.R. CASH, A BLOCK 6(4) RUNGE-KUTTA FORMULA FOR NON-STIFF C INITIAL VALUE PROBLEMS, C TO APPEAR IN TOMS. C C C C THE PARAMETERS REPRESENT- C FCN -- SUBROUTINE FCN(T,Y,YP) TO EVALUATE DERIVATIVES C YP(I)=DY(I)/DT C NEQN -- NUMBER OF DIFFERENTIAL EQUATIONS TO BE INTEGRATED. C Y(*) -- APPROXIMATION TO THE SOLUTION VECTOR AT T. C T -- INDEPENDENT VARIABLE. C TOUT -- THE NEXT POINT WHERE INTERMEDIATE OUTPUT IS C REQUIRED. APPROXIMATE SOLUTIONS AT THESE POINTS C WILL BE OBAINED BY INTERPOLATION. C TEND -- END OF THE INTEGRATION RANGE. THIS WILL BE HIT EXACTLY. C YP(*) -- APPROXIMATION TO THE DERIVATIVE VECTOR DY/DT AT T. C RELERR,ABSERR -- RELATIVE AND ABSOLUTE ERROR TOLERANCES FOR C LOCAL ERROR TEST. AT THE NTH POINT OF EACH BLOCK (N=1,2) C THE CODE REQUIRES THAT C DABS(LOCAL ERROR)/N .LE. RELERR*DABS(Y) + ABSERR C FOR EACH COMPONENT OF THE LOCAL ERROR AND SOLUTION VECTORS C IFLAG -- INDICATOR FOR STATUS OF INTEGRATION. C WORK(*) -- ARRAY TO HOLD INFORMATION INTERNAL TO BRKF45 WHICH C IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE DIMENSIONED C AT LEAST 6+9*NEQN C IWORK(*) -- INTEGER ARRAY USED TO HOLD INFORMATION INTERNAL TO C BRKF45 WHICH IS NECESSARY FOR SUBSEQUENT CALLS. MUST BE C DIMENSIONED AT LEAST 5 C C C*********************************************************************** C FIRST CALL TO BRKF45 C*********************************************************************** C C SUBROUTINE RKFC CONTAINS A MACHINE-DEPENDENT PARAMETER U. C THE SINGLE-PRECISION TOMS VERSION OF THIS PACKAGE CONTAINS A VALUE C SUITABLE FOR A CDC 6600, AND THE DOUBLE-PRECISION VERSION CONTAINS C A VALUE SUITABLE FOR AN IBM 3081. IF THE PACKAGE IS RUN ON ANY C OTHER MACHINE, THE USER SHOULD FIRST RESET U, WHICH IS DEFINED IN C A DATA STATEMENT, CLEARLY MARKED IN SUBROUTINE RKFC. APPROPRIATE C VALUES OF U FOR OTHER MAJOR MACHINES ARE SPECIFIED IN THE COMMENTS C OF RKFC. C C THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR THE ARRAYS C IN THE CALL LIST - Y(NEQN) , YP(NEQN) , WORK(6+9*NEQN) C IWORK(5), DECLARE FCN IN AN EXTERNAL STATEMENT, SUPPLY SUBROUTINE C FCN(T,Y,YP) AND INITIALIZE THE FOLLOWING PARAMETERS- C C NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED. (NEQN .GE. 1) C Y(*) -- VECTOR OF INITIAL CONDITIONS. C T -- STARTING POINT OF INTEGRATION , MUST BE A VARIABLE. C TOUT -- OUTPUT POINT AT WHICH SOLUTION, AND POSSIBLY THE C DERIVATIVE, IS DESIRED. C TEND -- END OF THE RANGE OF INTEGRATION. IF THE SOLUTION IS C REQUIRED ONLY AT TEND THEN THE USER SHOULD SET TOUT=TEND. C RELERR,ABSERR -- RELATIVE AND ABSOLUTE LOCAL ERROR TOLERANCES C WHICH MUST BE NON-NEGATIVE. RELERR MUST BE A VARIABLE WHILE C ABSERR MAY BE A CONSTANT. THE CODE SHOULD NORMALLY NOT BE C USED WITH RELATIVE ERROR CONTROL SMALLER THAN ABOUT 1.E-8, C UNLESS AN APPROPRIATE NONZERO ABSOLUTE TOLERANCE IS GIVEN. C TO AVOID LIMITING PRECISION DIFFICULTIES THE CODE REQUIRES C RELERR TO BE LARGER THAN AN INTERNALLY COMPUTED RELATIVE C ERROR PARAMETER WHICH IS MACHINE DEPENDENT. IN PARTICULAR, C PURE ABSOLUTE ERROR IS NOT PERMITTED. IF A SMALLER THAN C ALLOWABLE VALUE OF RELERR IS ATTEMPTED, BRKF45 INCREASES C RELERR APPROPRIATELY AND RETURNS CONTROL TO THE USER BEFORE C CONTINUING THE INTEGRATION. C IFLAG -- +1,-1 INDICATOR TO INITIALIZE THE CODE FOR EACH NEW C PROBLEM. NORMAL INPUT IS +1. THE USER SHOULD SET IFLAG=-1 C ONLY WHEN ONE-BLOCK INTEGRATOR CONTROL IS ESSENTIAL. IN C THIS CASE, BRKF45 ATTEMPTS TO ADVANCE THE SOLUTION A C SINGLE BLOCK IN THE DIRECTION OF TEND EACH TIME IT IS C CALLED. SINCE THIS MODE OF OPERATION RESULTS IN EXTRA C COMPUTING OVERHEAD, IT SHOULD BE AVOIDED UNLESS NEEDED. C C C*********************************************************************** C OUTPUT FROM BRKF45 C*********************************************************************** C C Y(*) -- COMPUTED SOLUTION APPROXIMATION AT T. C T -- VALUE OF THE INDEPENDENT VARIABLE WHERE THE SOLUTION IS C REPORTED. C IFLAG = 2 -- SUCCESSFUL RETURN. EITHER THE INTEGRATION REACHED C T=TEND OR A SUCCESSFUL INTERPOLATION HAS BEEN C PERFORMED AT T=TOUT. IF T.EQ.TEND THEN THE C INTEGRATION IS FINISHED. IF NOT, THE CODE SHOULD BE C CALLED WITH THE NEXT VALUE OF TOUT AND WITH IFLAG=+2 C FOR NORMAL INTEGRATION OR IFLAG=-2 FOR ONE-BLOCK C INTEGRATION. C =-2 -- A SINGLE SUCCESSFUL BLOCK IN THE DIRECTION OF TEND C HAS BEEN TAKEN. NORMAL MODE FOR CONTINUING C INTEGRATION ONE BLOCK AT A TIME. C = 3 -- INTEGRATION WAS NOT COMPLETED BECAUSE RELATIVE ERROR C TOLERANCE WAS TOO SMALL. RELERR HAS BEEN INCREASED C APPROPRIATELY FOR CONTINUING. C = 4 -- INTEGRATION WAS NOT COMPLETED BECAUSE MORE THAN C 18000 DERIVATIVE EVALUATIONS WERE NEEDED. THIS C IS APPROXIMATELY 2000 BLOCKS. C = 5 -- INTEGRATION WAS NOT COMPLETED BECAUSE SOLUTION C VANISHED MAKING A PURE RELATIVE ERROR TEST C IMPOSSIBLE. MUST USE NON-ZERO ABSERR TO CONTINUE. C USING THE ONE-BLOCK INTEGRATION MODE FOR ONE BLOCK C IS A GOOD WAY TO PROCEED. C = 6 -- INTEGRATION WAS NOT COMPLETED BECAUSE REQUESTED C ACCURACY COULD NOT BE ACHIEVED USING SMALLEST C ALLOWABLE STEPSIZE. USER MUST INCREASE THE ERROR C TOLERANCE BEFORE CONTINUED INTEGRATION CAN BE C ATTEMPTED. C = 7 -- INVALID INPUT PARAMETERS C THIS INDICATOR OCCURS IF ANY OF THE FOLLOWING IS C SATISFIED - NEQN .LE. 0 C T=TEND C RELERR OR ABSERR .LT. 0. C ABS(IFLAG) .LT. 1 OR .GT.7 C WORK(*),IWORK(*) -- INFORMATION WHICH IS USUALLY OF NO INTEREST C TO THE USER BUT NECESSARY FOR SUBSEQUENT CALLS. C WORK(1),...,WORK(NEQN) CONTAIN THE SOLUTION VECTOR C AND WORK(NEQN+1),...,WORK(2*NEQN) CONTAIN THE C DERIVATIVE VECTOR AT THE END POINT (WHICH IS ITSELF C CONTAINED IN WORK(2*NEQN+1)) OF THE BLOCK STEP JUST C COMPUTED. WORK(2*NEQN+2) CONTAINS THE STEPSIZE H C JUST USED. (THIS IS THE STEPSIZE BEING USED BY THE C INTERPOLANT OVER THIS BLOCK.) WORK(2*NEQN+3) C CONTAINS THE STEPSIZE H TO BE ATTEMPTED ON THE NEXT C BLOCK. IWORK(1) CONTAINS THE DERIVATIVE EVALUATION C COUNTER. C C C*********************************************************************** C SUBSEQUENT CALLS TO BRKF45 C*********************************************************************** C C SUBROUTINE BRKF45 RETURNS WITH ALL INFORMATION NEEDED TO CONTINUE C THE INTEGRATION. AFTER THE CODE REPORTS A SUCCESSFUL SOLUTION AT C TOUT (INDICATED BY IFLAG=2), THE USER NEEDS TO DEFINE A NEW TOUT C BEFORE SIMPLY CALLING BRKF45 AGAIN TO CONTINUE IN THE NORMAL MODE. C (BUT THE USER MUST FIRST RESET IFLAG TO -2 TO CONTINUE IN THE C ONE-BLOCK INTEGRATOR MODE.) C IF THE INTEGRATION WAS NOT COMPLETED BUT THE USER STILL WANTS TO C CONTINUE (IFLAG=3,4), HE JUST CALLS BRKF45 AGAIN. IN THE CASE C IFLAG=3 THE RELERR PARAMETER HAS BEEN ADJUSTED APPROPRIATELY FOR C CONTINUING THE INTEGRATION. IN THE CASE OF IFLAG=4 THE FUNCTION C COUNTER WILL BE RESET TO 0 AND ANOTHER 18000 FUNCTION EVALUATIONS C ARE ALLOWED. C HOWEVER,IN THE CASE IFLAG=5, THE USER MUST FIRST ALTER THE ERROR C CRITERION TO USE A POSITIVE VALUE OF ABSERR BEFORE INTEGRATION CAN C PROCEED. IF HE DOES NOT,EXECUTION IS TERMINATED. C ALSO,IN THE CASE IFLAG=6, IT IS NECESSARY FOR THE USER TO RESET C IFLAG TO 2 (OR -2 WHEN THE ONE-BLOCK INTEGRATION MODE IS BEING C USED) AS WELL AS INCREASING EITHER ABSERR,RELERR OR BOTH BEFORE C THE INTEGRATION CAN BE CONTINUED. IF THIS IS NOT DONE, EXECUTION C WILL BE TERMINATED. THE OCCURRENCE OF IFLAG=6 INDICATES A TROUBLE C SPOT(SOLUTION IS CHANGING RAPIDLY,SINGULARITY MAY BE PRESENT) AND C IT OFTEN IS INADVISABLE TO CONTINUE. C IF IFLAG=7 IS OBTAINED, INTEGRATION CAN NOT BE CONTINUED UNLESS C THE INVALID INPUT PARAMETERS ARE CORRECTED. C IT SHOULD BE NOTED THAT THE ARRAYS WORK,IWORK CONTAIN INFORMATION C REQUIRED FOR SUBSEQUENT INTEGRATION. ACCORDINGLY, WORK AND IWORK C SHOULD NOT BE ALTERED. C C C*********************************************************************** C USER CALLS TO THE INTERPOLANT ROUTINE EXTRA C*********************************************************************** C C SUBROUTINE EXTRA CAN ALSO BE CALLED BY THE USER, IN CONJUNCTION C WITH USAGE OF BRKF45, TO PROVIDE APPROXIMATE SOLUTIONS AT 'OFF-STEP C POINTS' BY USE OF THE INTERPOLATING POLYNOMIAL. WHILE BRKF45 C HANDLES THE USUAL SITUATIONS, IT CAN BE HELPFUL TO THE USER TO BE C ABLE TO ACCESS THE INTERPOLANT DIRECTLY, SUCH AS WHEN DOING ROOT C FINDING. ALSO, IT IS POSSIBLE THAT THE USER MAY HAVE SOME NEED FOR C EXTRAPOLATING OUTSIDE OF THE BLOCK STEP ON WHICH THE UNDERLYING C INTERPOLANT IS BASED. BRKF45 WILL NOT DO THIS. C THE FORM OF THE USAGE CALL IS C C CALL EXTRA ( NEQN, WORK(1), WORK(NEQN+1), WORK(2*NEQN+1), C WORK(2*NEQN+2), WORK(2*NEQN+4), WORK(3*NEQN+4), C WORK(4*NEQN+4), WORK(5*NEQN+4), TEX, YEX, YPEX ) C C WHERE YEX AND YPEX ARE THE SOLUTION VECTOR AND DERIVATIVE VECTOR C APPROXIMATIONS DEFINED BY THE INTERPOLANT AT THE POINT TEX, AND C WORK IS THE WORKING ARRAY SET UP BY BRKF45. C C C*********************************************************************** C MORE ABOUT THE INTERPOLANT ROUTINE EXTRA C*********************************************************************** C C THE FIRST DERIVATIVE APPROXIMATIONS GIVEN BY THE INTERPOLANT C CAN BE MADE LESS SENSITIVE TO ROUNDING ERRORS. IF VERY ACCURATE C FIRST DERIVATIVE APPROXIMATIONS ARE REQUIRED, THEN THREE C CHANGES SHOULD BE MADE TO THE CODE. THESE CHANGES ARE CLEARLY C DOCUMENTED IN THE CODE; JUST SEARCH FOR THE WORD 'CHANGES'. C (ONE LINE OF THE SUBROUTINE RKFC MUST BE REPLACED, ONE LINE C MUST BE ADDED TO THE SUBROUTINE RKFC, AND THE SUBROUTINE C EXTRA MUST BE REPLACED.) C FOR AN EXPLANATION OF THESE MODIFICATIONS SEE C 'REMARK ON ALGORITHM 669', BY D. HIGHAM; TO APPEAR IN TOMS. C C IF THE MODIFICATIONS ARE MADE, THEN THE CALL TO EXTRA DESCRIBED C ABOVE SHOULD BE CHANGED TO C C CALL EXTRA ( NEQN, WORK(6*NEQN+4), WORK(NEQN+1), WORK(2*NEQN+1), C WORK(2*NEQN+2), WORK(2*NEQN+4), WORK(3*NEQN+4), C WORK(4*NEQN+4), WORK(5*NEQN+4), TEX, YEX, YPEX ) C C C*********************************************************************** C IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL ENDPNT,BLKOUT DIMENSION Y(NEQN),YP(NEQN),WORK(*),IWORK(5) EXTERNAL FCN C C COMPUTE INDICES FOR THE SPLITTING OF THE WORK ARRAY C KW = 1 KWP = KW + NEQN KX = KWP + NEQN KHI = KX + 1 KH = KHI + 1 KY1 = KH + 1 KY2 = KY1 + NEQN KF1 = KY2 + NEQN KF2 = KF1 + NEQN KF3 = KF2 + NEQN KF4 = KF3 + NEQN KF7 = KF4 + NEQN KSR = KF7 + NEQN KSA = KSR + 1 KT = KSA + 1 C K ? = KT + 1 C THE WORK SPACE TOTALS 6 + 9*NEQN . C IF (IABS(IFLAG) .NE. 1) THEN ENDPNT = (IWORK(4) .EQ. -1) BLKOUT = (IWORK(5) .EQ. -1) END IF C C*********************************************************************** C THIS INTERFACING ROUTINE MERELY RELIEVES THE USER OF A LONG C CALLING LIST VIA THE SPLITTING APART OF TWO WORKING STORAGE C ARRAYS. IF THIS IS NOT COMPATIBLE WITH THE USERS COMPILER, C HE MUST USE RKFC DIRECTLY. C*********************************************************************** C CALL RKFC (FCN,NEQN,Y,T,TOUT,TEND,YP,RELERR,ABSERR,IFLAG,WORK(KW), * WORK(KWP),WORK(KX),WORK(KHI),WORK(KH),WORK(KY1),WORK(KY2), * WORK(KF1),WORK(KF2),WORK(KF3),WORK(KF4),WORK(KF7),WORK(KSR), * WORK(KSA),WORK(KT),IWORK(1),IWORK(2),IWORK(3),ENDPNT,BLKOUT) C IWORK(4) = 0 IF (ENDPNT) IWORK(4) = -1 IWORK(5) = 0 IF (BLKOUT) IWORK(5) = -1 C RETURN END SUBROUTINE RKFC (FCN,NEQN,Y,T,TOUT,TEND,YP,RELERR,ABSERR,IFLAG, * W,WP,X,HINT,H,Y1,Y2,F1,F2,F3,F4,F7,SAVRE,SAVAE, * TOLD,NFE,JFLAG,KFLAG,ENDPNT,BLKOUT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C TWO STEP BLOCK RUNGE-KUTTA FEHLBERG METHOD. C A STANDARD 5(4) FORMULA IS USED AT THE FIRST POINT IN THE BLOCK C AND A 6(4) FORMULA IS USED AT THE SECOND POINT. C LOGICAL HFAILD,ENDPNT,INTERP,BLKOUT DIMENSION Y(NEQN),YP(NEQN),Y1(NEQN),Y2(NEQN),F1(NEQN),F2(NEQN), * F3(NEQN),F4(NEQN),F7(NEQN),W(NEQN),WP(NEQN) EXTERNAL FCN C C*********************************************************************** C C COEFFICIENTS DEFINING THE METHOD ARE SET IN THE FOLLOWING C PARAMETER STATEMENTS. C PARAMETER( C2=0.25D+0, C3=3.0D+0/8.0D+0, C4=12.0D+0/13.0D+0, * C6=0.5D+0, C8=3.0D+0/2.0D+0, C9=2.0D+0 ) PARAMETER( A21=1.0D+0/4.0D+0) PARAMETER (A31=3.0D+0/32.0D+0, A32=9.0D+0/32.0D+0 ) PARAMETER( A41=1932.0D+0/2197.0D+0, A42=-7200.0D+0/2197.0D+0, * A43=7296.0D+0/2197.0D+0 ) PARAMETER( A51=439.0D+0/216.0D+0, A52=-8.0D+0, * A53=3680.0D+0/513.0D+0, A54=-845.0D+0/4104.0D+0) PARAMETER( A61=-8.0D+0/27.0D+0, A62=2.0D+0, * A63=-3544.0D+0/2565.0D+0, A64=1859.0D+0/4104.0D+0, * A65=-11.0D+0/40.0D+0) PARAMETER( B1B=16.0D+0/135.0D+0, B3B=6656.0D+0/12825.0D+0, * B4B=28561.0D+0/56430.0D+0, B5B=-9.0D+0/50.0D+0, * B6B=2.0D+0/55.0D+0) PARAMETER( ERC11=1.0D+0/360.0D+0, ERC13=-128.0D+0/4275.0D+0, * ERC14=-2197.0D+0/75240.0D+0, ERC15=1.0D+0/50.0D+0) C C THE ABOVE DEFINE THE COEFFICIENTS FOR BRKF45 USED TO GENERATE C THE SOLUTION AT THE FIRST BLOCK POINT. BELOW ARE ADDITIONAL C COEFFICIENTS NEEDED TO GENERATE THE SOLUTION AT THE SECOND BLOCK C POINT. C PARAMETER( B1=931.0D+0/6480.0D+0, B3=315392.0D+0/1500525.0D+0, * B4=371293.0D+0/615600.0D+0, B5=1.0D+0/50.0D+0, * B6=0.4D+0, B7=-4.0D+0/15.0D+0, * B8=85006.0D+0/115425.0D+0, B9=239.0D+0/1560.0D+0) PARAMETER( A81=-119397029895.0D+0/151948225000.0D+0, * A82=78390.0D+0/29081.0D+0, * A83=-51517464.0D+0/132821875.0D+0, * A84=-3780749193.0D+0/1168832500.0D+0, * A85=79268193.0D+0/55925000.0D+0, * A86=-11370591.0D+0/15379375.0D+0, * A87=5670.0D+0/2237.0D+0) PARAMETER( A91=23406188597.0D+0/8429231250.0D+0, * A92=-62928.0D+0/13623.0D+0, * A93=-31066887488.0D+0/5747203125.0D+0, * A94=164486461399.0D+0/8429231250.0D+0, * A95=-70336084.0D+0/11203125.0D+0, * A96=185680664.0D+0/24646875.0D+0, * A97=-3385330161.0D+0/243117160.0D+0, * A98=232648.0D+0/96795.0D+0) C C NEXT WE DEFINE COEFFICIENTS FOR THE ERROR ESTIMATE FORMULA. C PARAMETER(ERC21=0.22490537549467D-1,ERC23=-0.18660147916167D+0, * ERC24=0.224620349650850D+0, ERC25=-0.38462277720213D-1, * ERC26=9.0D+0/50.0D+0, ERC27=-7.0D+0/30.0D+0, * ERC28=0.36286203014893D-1, ERC29=-0.005D+0) C C IN WHAT FOLLOWS WE GIVE THE RATIONAL FORMS OF THESE DECIMAL C CONSTANTS THAT DEFINE THE ERROR COEFFICIENTS. C C B3E = 1067091077380.0D+0/1829119027671.0D+0 C B4E = 3284168845918.0D+0/21339721989495.0D+0 C B5E = 110317750789.0D+0/240996319200.0D+0- C * 4448925830089.0D+0/12329617149531.0D+0 C B6E = 1.0D+0/25.0D+0 C B7E = 0.2D+0 C B8E = 239992027043.0D+0/361494478800.0D+0 C B9E = 1273.0D+0/7800.0D+0 C B1E = 2.0D+0-B3E-B4E-B5E-B6E-B7E-B8E-B9E C ERC21 = (B1-B1E)/2.0D+0 C ERC23 = (B3-B3E)/2.0D+0 C ERC24 = (B4-B4E)/2.0D+0 C ERC25 = (B5-B5E)/2.0D+0 C ERC26 = 9.0D+0/50.0D+0 C ERC27 = -7.0D+0/30.0D+0 C ERC28 = (B8 - B8E)/2.0D+0 C ERC29 = -0.005D+0 C C IN WHAT FOLLOWS WE GIVE DECIMAL APPROXIMATIONS TO THE C COEFFICIENTS NEEDED TO COMPUTE THE SOLUTION OVER THE SECOND C STEP OF THE BLOCK. C C B1=0.14367283950617 C B3=0.21018776761467 C B4=0.60314002599090 C B5=0.02 C B6=0.4 C B7=-0.26666666666667 C B8=0.73646090534979 C B9=0.15320512820513 C A81=-0.78577443004023 C A82=2.6955744300402 C A83=-0.38786882055384 C A84=-3.2346372923408 C A85=1.4174017523469 C A86=-0.73934025277360 C A87=2.5346446133214 C A91=2.7767880489695 C A92=-4.6192468619247 C A93=-5.4055662923868 C A94=19.513815260318 C A95=-6.2782557545328 C A96=7.5336392138963 C A97=-13.924686192451 C A98=2.4035125781290 C B1E=0.098691764407238 C B3E=0.58339072593800 C B4E=0.1538993266892 C B5E=0.096924555440426 C B6E=0.04 C B7E=0.2 C B8E=0.66388849932001 C B9E=0.16320512820513 C C THE COEFFICIENTS OF THE RUNGE-KUTTA FORMULA ARE NOW SET. C C*********************************************************************** C C THE COMPUTER UNIT ROUNDOFF ERROR U IS THE SMALLEST POSITIVE VALUE C REPRESENTABLE IN THE MACHINE SUCH THAT 1.+ U .GT. 1. C VALUES TO BE USED ARE C U = 9.5E-7 FOR IBM 360/370 (SINGLE PRECISION) C U = 1.2E-7 FOR DEC VAX (SINGLE PRECISION) C U = 1.5E-8 FOR UNIVAC 1108 C U = 7.5E-9 FOR PDP-10 C U = 7.1D-15 FOR CDC 6000 SERIES C U = 2.2D-16 FOR IBM 360/370 (DOUBLE PRECISION) C U = 2.8D-17 FOR DEC VAX (DOUBLE PRECISION) C C DATA U / 1.2D-7 / DATA U / 2.2D-16 / C C*********************************************************************** C C REMIN IS A TOLERANCE THRESHOLD WHICH IS ALSO DETERMINED BY THE C INTEGRATION METHOD. IN PARTICULAR, A FIFTH ORDER METHOD WILL C GENERALLY NOT BE CAPABLE OF DELIVERING ACCURACIES NEAR LIMITING C PRECISION ON COMPUTERS WITH LONG WORDLENGTHS. THIS DOES NOT HAVE TO C BE CHANGED FOR DIFFERENT MACHINES. C DATA REMIN / 1.0D-13 / C C*********************************************************************** C C THE EXPENSE IS CONTROLLED BY RESTRICTING THE NUMBER C OF FUNCTION EVALUATIONS TO BE APPROXIMATELY MAXNFE. C DATA MAXNFE / 18000 / C C*********************************************************************** C C C CHECK INPUT PARAMETERS C C IF (NEQN.LT.1) GO TO 10 IF ((RELERR.LT.0.0D+0).OR.(ABSERR.LT.0.0D+0)) GO TO 10 IF (T.EQ.TEND) GO TO 10 MFLAG = IABS(IFLAG) IF (MFLAG .EQ. 1) THEN TOLD = T IF (TOUT .EQ. T) GO TO 20 IF (DSIGN(1.0D+0,TEND-T) .NE. DSIGN(1.0D+0,TOUT-T)) GO TO 10 ELSE IF (TOUT .EQ. TOLD) GO TO 20 IF (DSIGN(1.0D+0,TEND-TOLD) .NE. DSIGN(1.0D+0,TOUT-TOLD)) * GO TO 10 END IF IF ((MFLAG.GE.1).AND.(MFLAG.LE.7)) GO TO 20 C C INVALID INPUT, RETURN C 10 IFLAG = 7 RETURN C C IS THIS THE FIRST CALL? IF SO JUMP TO 70. C 20 IF (MFLAG.EQ.1) GO TO 70 C C CHECK CONTINUATION POSSIBILITIES C IF (MFLAG.NE.2) GO TO 30 C C IFLAG = +2 OR -2 C IF (KFLAG.LT.3) GO TO 70 IF (KFLAG.EQ.3) GO TO 60 IF (KFLAG.EQ.4) GO TO 50 IF ((KFLAG.EQ.5).AND.(ABSERR.EQ.0.0D+0)) GO TO 40 IF ((KFLAG.EQ.6).AND.(RELERR.LE.SAVRE).AND.(ABSERR.LE.SAVAE)) * GO TO 40 GO TO 70 C C IFLAG = 3,4,5,6 OR 7 C 30 IF (IFLAG.EQ.3) GO TO 60 IF (IFLAG.EQ.4) GO TO 50 IF ((IFLAG.EQ.5).AND.(ABSERR.GT.0.0D+0)) GO TO 60 C C INTEGRATION CANNOT BE CONTINUED SINCE USER DID NOT RESPOND TO C THE INSTRUCTIONS PERTAINING TO IFLAG=5,6 OR 7 C 40 STOP C C C*********************************************************************** C C RESET FUNCTION EVALUATION COUNTER C 50 NFE = 0 IF (MFLAG.EQ.2) GO TO 70 C C RESET FLAG VALUE FROM PREVIOUS CALL C 60 IFLAG = JFLAG IF (KFLAG.EQ.3) MFLAG = IABS(IFLAG) C C SAVE INPUT IFLAG AND SET CONTINUATION FLAG VALUE FOR SUBSEQUENT C INPUT CHECKING C 70 JFLAG = IFLAG KFLAG = 0 C C SAVE RELERR AND ABSERR FOR CHECKING INPUT ON SUBSEQUENT CALLS C SAVRE = RELERR SAVAE = ABSERR C C RESTRICT RELATIVE ERROR TOLERANCE TO BE AT LEAST AS LARGE AS C 2U+REMIN TO AVOID LIMITING PRECISION DIFFICULTIES ARISING FROM C IMPOSSIBLE ACCURACY REQUESTS. IF TOLERANCE TOO SMALL, INCREASE C AND RETURN. C RER = 2.0D+0*U+REMIN IF (RELERR .LT. RER) THEN C C RELATIVE ERROR TOLERANCE TOO SMALL C RELERR = RER IFLAG = 3 KFLAG = 3 RETURN END IF C U26 = 26.0D+0*U IF (MFLAG.NE.1) GO TO 100 C C*********************************************************************** C C INITIALIZATION -- C DEFINE INTEGRATION INDEPENDENT VARIABLE X C EVALUATE INITIAL DERIVATIVES C SET UP WORKING ARRAYS FOR INTEGRATION VARIABLES C SET COUNTER FOR FUNCTION EVALUATIONS,NFE C ESTIMATE STARTING STEPSIZE C X = T ENDPNT = .FALSE. BLKOUT = .FALSE. C A = T CALL FCN (A,Y,YP) NFE = 1 DO 80 N = 1, NEQN W(N) = Y(N) WP(N) = YP(N) 80 CONTINUE C C COMPUTE INITIAL STEPLENGTH. C DT = TOUT - T IF (DT .EQ. 0.0D+0) DT = TEND - T H = DABS(DT) TOLN = 0.0D+0 DO 90 K = 1, NEQN TOL = RELERR*DABS(Y(K))+ABSERR IF (TOL.LE.0.0D+0) GO TO 90 TOLN = TOL YPK = DABS(YP(K)) IF (YPK*H**5.GT.TOL) H = (TOL/YPK)**0.2D+0 90 CONTINUE IF (TOLN.LE.0.0D+0) H = 0.0D+0 H = DMAX1(H,U26*DMAX1(DABS(T),DABS(DT))) JFLAG = ISIGN(2,IFLAG) H = DSIGN(H,DT) C C INITIAL STEPLENGTH NOW COMPUTED. COMPUTE FIRST SOLUTION. C C*********************************************************************** C 100 CONTINUE C C TO AVOID PREMATURE UNDERFLOW IN THE ERROR TOLERANCE FUNCTION, C SCALE THE ERROR TOLERANCES. C SCALE = 2.0D+0/RELERR AE = SCALE*ABSERR C C SET SAFETY FACTOR FOR STEPSIZE ADJUSTMENT, BASED ON TOLERANCES. C TOLER = DMAX1(ABSERR,RELERR) SF = 0.85D+0 IF (TOLER .GE. 1.D-5) SF = 0.8D+0 IF (TOLER .LE. 1.D-9) SF = 0.9D+0 C C RESTORE INTEGRATION VARIABLE TO END OF LAST BLOCK STEP TAKEN C AND SET THE SIGN OF THE DIRECTION OF INTEGRATION. C T = X DTSIGN = DSIGN(1.0D+0,TEND-T) C C HAVE WE ALREADY INTEGRATED PAST THE PRESENT DATA OUTPUT POINT? C IF SO JUMP TO 390 AND PERFORM INTERPOLATION. C IF NOT, SEE IF WE HAVE REACHED THE END POINT OF INTEGRATION. C IF NOT, SEE IF RESULTS AT THE END OF THE BLOCK STEP NEED TO BE C REPORTED. C IF ((TOUT-T)*DTSIGN .LE. 0.0D+0) THEN IF (TOUT .EQ. T) THEN IFLAG = 2 ENDPNT = .FALSE. GO TO 400 ELSE GO TO 390 END IF END IF IF (ENDPNT) THEN IFLAG = 2 ENDPNT = .FALSE. GO TO 400 END IF IF (IFLAG .EQ. -2 .AND. BLKOUT) THEN BLKOUT = .FALSE. GO TO 400 END IF C C*********************************************************************** C*********************************************************************** C BLOCK BY BLOCK INTEGRATION C 150 CONTINUE C C SEE IF WE ARE TOO CLOSE TO THE END POINT. IF SO, DO LINEAR C EXTRAPOLATION AND RETURN. C IF (DABS(TEND-T) .LE. U26*DABS(T)) THEN IF (DTSIGN*(TEND-TOUT) .GT. 0.0D+0) THEN DT = TOUT - T ENDPNT = .FALSE. T = TOUT ELSE DT = TEND - T ENDPNT = .TRUE. T = TEND X = TEND END IF DO 160 K = 1, NEQN Y(K) = W(K)+DT*WP(K) 160 CONTINUE CALL FCN (T,Y,YP) NFE = NFE+1 IF (ENDPNT) THEN DO 170 N=1,NEQN W(N) = Y(N) WP(N) = YP(N) 170 CONTINUE ENDPNT = .FALSE. END IF IFLAG = 2 RETURN END IF C C C SET SMALLEST ALLOWABLE STEPSIZE AND STEP FAILURE FLAG. C ADJUST STEPSIZE IF NECESSARY TO HIT THE END POINT OF INTEGRATION. C HMIN = U26*DABS(T) HFAILD = .FALSE. HSTOP = 0.5D+0*(TEND-T) IF (DABS(HSTOP) .LE. DABS(H)) THEN ENDPNT = .TRUE. H = HSTOP END IF C C*********************************************************************** C CORE INTEGRATOR FOR A SINGLE BLOCK C*********************************************************************** C THE TOLERANCES HAVE BEEN SCALED TO AVOID PREMATURE UNDERFLOW IN C COMPUTING THE ERROR TOLERANCE FUNCTION ERRTOL. C TO AVOID PROBLEMS WITH ZERO CROSSINGS,RELATIVE ERROR IS MEASURED C USING THE AVERAGE OF THE MAGNITUDES OF THE SOLUTION AT THE C BEGINNING AND END POINTS OF A BLOCK. C TO DISTINGUISH THE VARIOUS ARGUMENTS, H IS NOT PERMITTED C TO BECOME SMALLER THAN 26 UNITS OF ROUNDOFF IN T. C PRACTICAL LIMITS ON THE CHANGE IN THE STEPSIZE ARE ENFORCED TO C SMOOTH THE STEPSIZE SELECTION PROCESS AND TO AVOID EXCESSIVE C CHATTERING ON PROBLEMS HAVING DISCONTINUITIES. C*********************************************************************** C C C TEST NUMBER OF DERIVATIVE FUNCTION EVALUATIONS. C IF OKAY,TRY TO ADVANCE THE INTEGRATION FROM T TO T+H C 200 IF (NFE .GT. MAXNFE) THEN C C TOO MUCH WORK C IFLAG = 4 KFLAG = 4 GO TO 400 END IF C C ADVANCE AN APPROXIMATE SOLUTION OVER ONE STEP OF LENGTH H C DO 210 N = 1, NEQN Y1(N) = W(N)+A21*H*WP(N) 210 CONTINUE CALL FCN (T+C2*H,Y1,F2) DO 220 N = 1, NEQN Y1(N) = W(N)+H*(A31*WP(N)+A32*F2(N)) 220 CONTINUE CALL FCN (T+C3*H,Y1,F3) DO 230 N = 1, NEQN Y1(N) = W(N)+H*(A41*WP(N)+A42*F2(N)+A43*F3(N)) 230 CONTINUE CALL FCN (T+C4*H,Y1,F4) DO 240 N = 1, NEQN Y1(N) = W(N)+H*(A51*WP(N)+A52*F2(N)+A53*F3(N)+A54*F4(N)) 240 CONTINUE CALL FCN (T+H,Y1,Y) DO 250 N = 1, NEQN Y1(N) = W(N)+H*(A61*WP(N)+A62*F2(N)+A63*F3(N)+A64*F4(N)+ * A65*Y(N)) 250 CONTINUE CALL FCN (T+C6*H,Y1,YP) NFE = NFE+5 EEOET = 0.0D+0 DO 270 N = 1, NEQN Y1(N) = W(N)+H*(B1B*WP(N)+B3B*F3(N)+B4B*F4(N)+B5B*Y(N)+ * B6B*YP(N)) ERRTOL = DABS(W(N))+DABS(Y1(N))+AE IF (ERRTOL .LE. 0.0D+0) THEN C C INAPPROPRIATE ERROR TOLERANCE C IFLAG = 5 KFLAG = 5 GO TO 400 END IF EZ = DABS(H*(ERC11*WP(N)+ERC13*F3(N)+ERC14*F4(N)+ERC15*Y(N)+ * B6B*YP(N))) EEOET = DMAX1(EEOET,EZ/ERRTOL) 270 CONTINUE ESTTOL = EEOET*SCALE C C CHECK THE ERROR ESTIMATE. IF STEPLENGTH HAS FAILED FOR FIRST STEP C IN THE BLOCK, GO AND COMPUTE A SMALLER STEP FOR A RE-TRY. C IF (ESTTOL .GT. 1.0D+0) GO TO 320 C C C INTEGRATION SUCCESSFUL FOR FIRST STEP. NOW INTEGRATE OVER C THE SECOND STEP IN THE BLOCK C CALL FCN (T+H,Y1,F1) DO 290 N = 1, NEQN Y2(N) = W(N)+H*(A81*WP(N)+A82*F2(N)+A83*F3(N)+A84*F4(N)+ * A85*Y(N)+A86*YP(N)+A87*F1(N)) 290 CONTINUE CALL FCN (T+C8*H,Y2,F7) DO 300 N = 1, NEQN Y2(N) = W(N)+H*(A91*WP(N)+A92*F2(N)+A93*F3(N)+A94*F4(N)+ * A95*Y(N)+A96*YP(N)+A97*F1(N)+A98*F7(N)) 300 CONTINUE CALL FCN (T+C9*H,Y2,F2) NFE = NFE+3 EEOET = 0.0D+0 DO 310 N = 1, NEQN EE = DABS(H*(ERC21*WP(N)+ERC23*F3(N)+ERC24*F4(N)+ERC25*Y(N)+ * ERC26*YP(N)+ERC27*F1(N)+ERC28*F7(N)+ERC29*F2(N))) Y2(N) = W(N)+H*(B1*WP(N)+B3*F3(N)+B4*F4(N)+B5*Y(N)+B6*YP(N)+ * B7*F1(N)+B8*F7(N)+B9*F2(N)) ERRTOL = DABS(Y2(N))+DABS(Y1(N))+AE IF (ERRTOL .LE. 0.0D+0) THEN C C INAPPROPRIATE ERROR TOLERANCE C IFLAG = 5 KFLAG = 5 GO TO 400 END IF EEOET = DMAX1(EEOET,EE/ERRTOL) 310 CONTINUE ESTTOL = EEOET*SCALE C C THE FIRST OF THE THREE CHANGES NEEDED TO MAKE THE INTERPOLANT C DERIVATIVE MORE STABLE IS TO UNCOMMENT THE FOLLOWING DO-LOOP C C DO 315 N = 1, NEQN C Y1(N) = B1B*WP(N)+B3B*F3(N)+B4B*F4(N)+B5B*Y(N)+ C * B6B*YP(N) C F3(N) = B1*WP(N)+B3*F3(N)+B4*F4(N)+B5*Y(N)+B6*YP(N)+ C * B7*F1(N)+B8*F7(N)+B9*F2(N) C 315 CONTINUE C C CHECK THE ERROR ESTIMATE OVER THE SECOND STEP OF THE BLOCK. C IF (ESTTOL .LE. 1.0D+0) GO TO 330 C C UNSUCCESSFUL BLOCK C REDUCE THE STEPSIZE , TRY AGAIN C THE DECREASE IS LIMITED TO A FACTOR OF ABOUT 1/10 C 320 HFAILD = .TRUE. S = 0.1D+0 ENDPNT = .FALSE. IF (ESTTOL.LT.1.0D+5) S = SF/ESTTOL**0.2D+0 H = S*H IF (DABS(H) .LT. HMIN) THEN C C REQUESTED ERROR UNATTAINABLE AT SMALLEST ALLOWABLE STEPSIZE C IFLAG = 6 KFLAG = 6 GO TO 400 END IF GO TO 200 C C C SUCCESSFUL BLOCK. CHECK FOR NEED TO INTERPOLATE. C STORE SOLUTION AT T+2*H. C 330 T = T+2.0D+0*H IF (ENDPNT) T = TEND TOLD = X X = T HINT = H INTERP = .FALSE. IF ((T-TOUT)*DTSIGN .GE. 0.0D+0) INTERP = .TRUE. IF (INTERP) THEN DO 340 N=1,NEQN SWAP = W(N) W(N) = Y2(N) Y2(N) = SWAP F2(N) = WP(N) 340 CONTINUE ELSE DO 350 N=1,NEQN W(N) = Y2(N) 350 CONTINUE END IF C A = T CALL FCN (A,W,WP) NFE = NFE+1 C C CHOOSE NEXT STEPSIZE C THE INCREASE IS LIMITED TO A FACTOR OF ABOUT 10 C IF STEP FAILURE HAS JUST OCCURRED, NEXT STEP IS NOT ALLOWED C TO INCREASE. C S = 10.0D+0 IF (ESTTOL.GT.1.0D-5) S = SF/ESTTOL**0.2D+0 IF (HFAILD) S = DMIN1(S,1.0D+0) H = DSIGN(DMAX1(S*DABS(H),HMIN),H) C C HAVE WE INTEGRATED PAST AN OUTPUT POINT? C IF SO, CALL THE INTERPOLATION ROUTINE AT T=TOUT. C IF NOT, SEE IF WE ARE AT THE END POINT OF INTEGRATION. C OTHERWISE, CHECK IF USER WANTS SOLUTIONS AT THE END OF THE BLOCK C STEP, OR ELSE CONTINUE THE INTEGRATION. C IF (INTERP) GO TO 390 IF (ENDPNT) THEN IFLAG = 2 ENDPNT = .FALSE. GO TO 400 END IF IF (IFLAG .LT. 0) THEN IFLAG = -2 GO TO 400 ELSE GO TO 150 END IF C C INTERPOLATE TO GET DATA AT OFF-STEP POINT AND RETURN. C C THE SECOND OF THE THREE CHANGES NEEDED TO MAKE THE INTERPOLANT C DERIVATIVE MORE STABLE IS TO REPLACE STATEMENT 390 BY THE C VERSION IN THE COMMENT STATEMENT BELOW C 390 CALL EXTRA (NEQN,W,WP,X,HINT,Y1,Y2,F1,F2,TOUT,Y,YP) C 390 CALL EXTRA (NEQN,F3,WP,X,HINT,Y1,Y2,F1,F2,TOUT,Y,YP) T = TOUT IF (IFLAG .LT. 0 .AND. TOUT .NE. X) BLKOUT = .TRUE. IFLAG = 2 RETURN C C RETURN WITH THE SOLUTION AT THE END OF THE LAST SUCCESSFUL C BLOCK STEP. C 400 DO 410 N=1,NEQN Y(N) = W(N) YP(N) = WP(N) 410 CONTINUE RETURN C END C C THE THIRD OF THE THREE CHANGES NEEDED TO MAKE THE INTERPOLANT C DERIVATIVE MORE STABLE IS TO REPLACE THE SUBROUTINE EXTRA BY C THE VERSION IN COMMENT STATEMENTS BELOW C SUBROUTINE EXTRA (NEQN,Y2,F2,T,H,Y1,Y,F1,F,TEX,YEX,YPEX) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C T IS THE INDEPENDENT VARIABLE; H IS THE STEP SPACING FOR VALUES C Y,Y1,Y2; TEX IS THE POINT WHERE THE OUTPUT IS REQUIRED; C F,F1 AND F2 ARE THE DERIVATIVES AT T-2H, T-H AND T. C YEX WILL HOLD THE APPROXIMATE SOLUTION AT TEX AND YPEX WILL C HOLD THE DERIVATIVE APPROXIMATION AT THIS POINT. C DIMENSION Y(NEQN),Y1(NEQN),Y2(NEQN),F(NEQN),F1(NEQN),F2(NEQN), * YEX(NEQN),YPEX(NEQN) C C PERFORM QUINTIC INTERPOLATION BASED ON THE VALUES C Y(N),Y1(N),Y2(N),F(N),F1(N),F2(N) AT THE POINTS C T-2H, T-H, T. THIS POLYNOMIAL IS EVALUATED AT THE POINT TEX C TO OBTAIN THE APPROXIMATE VALUE OF Y(TEX) AND THIS IS STORED C IN THE ARRAY YEX. THE DERIVATIVE OF THE INTERPOLATING POLYNOMIAL C IS ALSO EVALUATED AT TEX TO OBTAIN AN APPROXIMATION TO DY/DT C AT THIS POINT AND THIS APPROXIMATION IS STORED IN YPEX. C C SIG = (TEX-T)/H+2.0D+0 DO 10 N = 1, NEQN HF = H*F(N) HF1 = H*F1(N) HF2 = H*F2(N) A1 = 2.0D+0*(HF2+HF)+6.0D+0*(Y(N)-Y2(N))+8.0D+0*HF1 A2 = Y2(N)+4.0D+0*(Y1(N)-HF1)-5.0D+0*Y(N)-2.0D+0*HF A3 = HF1+HF+2.0D+0*(Y(N)-Y1(N)) A4 = Y1(N)-Y(N)-HF YEX(N) = ((((A1*0.125D+0*(SIG-2.0D+0)+0.25D+0*A2)*(SIG-1.0D+0)+ * A3)*(SIG-1.0D+0)+A4)*SIG+HF)*SIG+Y(N) YPEX(N) = ((((0.625D+0*A1*SIG+(A2-2.0D+0*A1))*SIG+(1.875D+0*A1- * 1.5D+0*A2+3.0D+0*A3))*SIG+(0.5D+0*(A2-A1)- * 2.0D+0*(A3-A4)))*SIG+HF)/H 10 CONTINUE RETURN END C C C C SUBROUTINE EXTRA (NEQN,YINC2,F2,T,H,YINC1,Y,F1,F,TEX,YEX,YPEX) C IMPLICIT DOUBLE PRECISION(A-H,O-Z) CC CC T IS THE INDEPENDENT VARIABLE; CC Y IS THE VALUE AT T-2H, THE VALUES AT T-H AND T ARE CC Y + H*YINC1 AND Y + H*YINC2 RESPECTIVELY; CC F,F1 AND F2 ARE THE DERIVATIVES AT T-2H, T-H AND T. CC TEX IS THE POINT WHERE THE OUTPUT IS REQUIRED; CC YEX WILL HOLD THE APPROXIMATE SOLUTION AT TEX AND YPEX WILL CC HOLD THE DERIVATIVE APPROXIMATION AT THIS POINT. CC C DIMENSION Y(NEQN),YINC1(NEQN),YINC2(NEQN),F(NEQN),F1(NEQN), C * F2(NEQN),YEX(NEQN),YPEX(NEQN) CC CC PERFORM QUINTIC INTERPOLATION BASED ON THE VALUES CC Y(N),Y(N)+H*YINC1(N),Y(N)+H*YINC2(N),F(N),F1(N),F2(N) AT THE CC POINTS T-2H, T-H, T. THIS POLYNOMIAL IS EVALUATED AT THE POINT CC TEX TO OBTAIN THE APPROXIMATE VALUE OF Y(TEX) AND THIS IS STORED CC IN THE ARRAY YEX. THE DERIVATIVE OF THE INTERPOLATING POLYNOMIAL CC IS ALSO EVALUATED AT TEX TO OBTAIN AN APPROXIMATION TO DY/DT CC AT THIS POINT AND THIS APPROXIMATION IS STORED IN YPEX. CC CC C SIG = (TEX-T)/H+2.0D+0 C DO 10 N = 1, NEQN C HF = H*F(N) C HF1 = H*F1(N) C HF2 = H*F2(N) C A1 = 2.0D+0*(HF2+HF)-6.0D+0*H*YINC2(N)+8.0D+0*HF1 C A2 = H*YINC2(N)+4.0D+0*H*YINC1(N)-4.0D+0*HF1-2.0D+0*HF C A3 = HF1+HF-2.0D+0*H*YINC1(N) C A4 = H*YINC1(N)-HF C YEX(N) = ((((A1*0.125D+0*(SIG-2.0D+0)+0.25D+0*A2)*(SIG-1.0D+0)+ C * A3)*(SIG-1.0D+0)+A4)*SIG+HF)*SIG+Y(N) C YPEX(N) = ((((0.625D+0*A1*SIG+(A2-2.0D+0*A1))*SIG+(1.875D+0*A1- C * 1.5D+0*A2+3.0D+0*A3))*SIG+(0.5D+0*(A2-A1)- C * 2.0D+0*(A3-A4)))*SIG+HF)/H C 10 CONTINUE C RETURN C END