C ALGORITHM 407 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN COMM. ACM, VOL. 14, NO. 03, C P. 185. SUBROUTINE DIFSUB(N,T,Y,SAVE,H,HMIN,HMAX,EPS,MF,YMAX,ERROR,KFLAG, 1 JSTART,MAXDER,PW) DOUBLE PRECISION A,D,E,H,R,T,Y,R1,R2,BND,EPS,EUP,EDWN,ENQ1 1 ,ENQ2,ENQ3,HMAX,HMIN,HNEW,HOLD,SAVE,TOLD,YMAX,ERROR,RACUM C 50,80,1010,1160,1230,1380,1410,1560,1580,1640,1710,1790,1820, C 1900,1970,2730,2770,2840,2920,2970,3020,3560,3630,3780,3950, C 4870,4900 ARE BLANK COMMENT CARDS IN THE PUBLISHED ALGORITHM C ACCORDING TO PUBLISHED ALGORITHM DFLOAT(K) IS CORRECT, NOT C DBLE(FLOAT(K)) ON LINES 2700,4390 C BUT DFLOAT(K) IS FOUND AS UNSATISFIED EXTERNAL REFERENCE WHILE C COMPILING ON CDC C THE PARAMETERS TO THE SUBROUTINE DIFSUB HAVE C THE FOLLOWING MEANINGS.. C N THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. N C MAY BE DECREASED ON LATER CALLS IF THE NUMBER OF C ACTIVE EQUATIONS REDUCES, BUT IT MUST NOT BE C INCREASED WITHOUT CALLING WITH JSTART = 0. C T THE INDEPENDENT VARIABLE. C Y AN 8 BY N ARRAY CONTAINING THE DEPENDENT VARIABLES AND C THEIR SCALED DERIVATIVES. Y(J+1,I) CONTAINS C THE J-TH DERIVATIVE OF Y(I) SCALED BY C H**J/FACTORIAL(J) WHERE H IS THE CURRENT C STEP SIZE. ONLY Y(1,I) NEED BE PROVIDED BY C THE CALLING PROGRAM ON THE FIRST ENTRY. C IF IT IS DESIRED TO INTERPOLATE TO NON MESH POINTS C THESE VALUES CAN BE USED. IF THE CURRENT STEP SIZE C IS H AND THE VALUE AT T + E IS NEEDED, FORM C S = E/H, AND THEN COMPUTE C NQ C Y(I)(T+E) = SUM Y(J+1,I)*S**J C J=0 C SAVE A BLOCK OF AT LEAST 12*N FLOATING POINT LOCATIONS C USED BY THE SUBROUTINES. C H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. C H MAY BE ADJUSTED UP OR DOWN BY THE PROGRAM C IN ORDER TO ACHIEVE AN ECONOMICAL INTEGRATION. C HOWEVER, IF THE H PROVIDED BY THE USER DOES C NOT CAUSE A LARGER ERROR THAN REQUESTED, IT C WILL BE USED. TO SAVE COMPUTER TIME, THE USER IS C ADVISED TO USE A FAIRLY SMALL STEP FOR THE FIRST C CALL. IT WILL BE AUTOMATICALLY INCREASED LATER. C HMIN THE MINIMUM STEP SIZE THAT WILL BE USED FOR THE C INTEGRATION. NOTE THAT ON STARTING THIS MUST C MUCH SMALLER THAN THE AVERAGE H EXPECTED SINCE C A FIRST ORDER METHOD IS USED INITIALLY. C HMAX THE MAXIMUM SIZE TO WHICH THE STEP WILL BE INCREASED C EPS THE ERROR TEST CONSTANT. SINGLE STEP ERROR ESTIMATES C DIVIDED BY YMAX(I) MUST BE LESS THAN THIS C IN THE EUCLIDEAN NORM. THE STEP AND/OR ORDER IS C ADJUSTED TO ACHIEVE THIS. C MF THE METHOD INDICATOR. THE FOLLOWING ARE ALLOWED.. C 0 AN ADAMS PREDICTOR CORRECTOR IS USED. C 1 A MULTI-STEP METHOD SUITABLE FOR STIFF C SYSTEMS IS USED. IT WILL ALSO WORK FOR C NON STIFF SYSTEMS. HOWEVER THE USER C MUST PROVIDE A SUBROUTINE PEDERV WHICH C EVALUATES THE PARTIAL DERIVATIVES OF C THE DIFFERENTIAL EQUATIONS WITH RESPECT C TO THE Y*S. THIS IS DONE BY CALL C PEDERV(T,Y,PW,M). PW IS AN N BY N ARRAY C WHICH MUST BE SET TO THE PARTIAL OF C THE I-TH EQUATION WITH RESPECT C TO THE J DEPENDENT VARIABLE IN PW(I,J). C PW IS ACTUALLY STORED IN AN M BY M C ARRAY WHERE M IS THE VALUE OF N USED ON C THE FIRST CALL TO THIS PROGRAM. C 2 THE SAME AS CASE 1, EXCEPT THAT THIS C SUBROUTINE COMPUTES THE PARTIAL C DERIVATIVES BY NUMERICAL DIFFERENCING C OF THE DERIVATIVES. HENCE PEDERV IS C NOT CALLED. C YMAX AN ARRAY OF N LOCATIONS WHICH CONTAINS THE MAXIMUM C OF EACH Y SEEN SO FAR. IT SHOULD NORMALLY BE SET TO C 1 IN EACH COMPONENT BEFORE THE FIRST ENTRY. (SEE THE C DESCRIPTION OF EPS.) C ERROR AN ARRAY OF N ELEMENTS WHICH CONTAINS THE ESTIMATED C ONE STEP ERROR IN EACH COMPONENT. C KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. C +1 THE STEP WAS SUCCESFUL. C -1 THE STEP WAS TAKEN WITH H = HMIN, BUT THE C REQUESTED ERROR WAS NOT ACHIEVED. C -2 THE MAXIMUM ORDER SPECIFIED WAS FOUND TO C BE TOO LARGE. C -3 CORRECTOR CONVERGENCE COULD NOT BE C ACHIEVED FOR H .GT. HMIN. C -4 THE REQUESTED ERROR IS SMALLER THAN CAN C BE HANDLED FOR THIS PROBLEM. C JSTART AN INPUT INDICATOR WITH THE FOLLOWING MEANINGS.. C -1 REPEAT THE LAST STEP WITH A NEW H C 0 PERFORM THE FIRST STEP. THE FIRST STEP C MUST BE DONE WITH THIS VALUE OF JSTART C SO THAT THE SUBROUTINE CAN INITIALIZE C ITSELF. C +1 TAKE A NEW STEP CONTINUING FROM THE LAST. C JSTART IS SET TO NQ, THE CURRENT ORDER OF THE METHOD C AT EXIT. NQ IS ALSO THE ORDER OF THE MAXIMUM C DERIVATIVE AVAILABLE. C MAXDER THE MAXIMUM DERIVATIVE THAT SHOULD BE USED IN THE C METHOD. SINCE THE ORDER IS EQUAL TO THE HIGHEST C DERIVATIVE USED, THIS RESTRICTS THE ORDER. IT MUST C BE LESS THAN 8 OR 7 FOR ADAMS OR STIFF METHODS C RESPECTIVELY. C PW A BLOCK OF AT LEAST N**2 FLOATING POINT LOCATIONS. DIMENSION Y(8,N),YMAX(N),SAVE(10,N),ERROR(N),PW(N), 1 A(8),PERTST(7,2,3) C THE COEFFICIENTS IN PERTST ARE USED IN SELECTING THE STEP AND C ORDER, THEREFORE ONLY ABOUT ONE PERCENT ACCURACY IS NEEDED. DATA PERTST /2.0,4.5,7.333,10.42,13.7,17.15,1.0, 1 2.0,12.0,24.0,37.89,53.33,70.08,87.97, 1 3.0,6.0,9.167,12.5,15.98,1.0,1.0, 1 12.0,24.0,37.89,53.33,70.08,87.97,1.0, 1 1.,1.,0.5,0.1667,0.04133,0.008267,1.0, 1 1.0,1.0,2.0,1.0,.3157,.07407,.0139/ DATA A(2) / -1.0/ IRET = 1 KFLAG = 1 IF (JSTART.LE.0) GO TO 140 C BEGIN BY SAVING INFORMATION FOR POSSIBLE RESTARTS AND CHANGING C H BY THE FACTOR R IF THE CALLER HAS CHANGED H. ALL VARIABLES C DEPENDENT ON H MUST ALSO BE CHANGED. C E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS C TO TEST FOR INCREASING THE ORDER, EDWN FOR DECREASING THE ORDER. C HNEW IS THE STEP SIZE THAT WAS USED ON THE LAST CALL. 100 DO 110 I = 1,N DO 110 J = 1,K 110 SAVE(J,I) = Y(J,I) HOLD = HNEW IF (H.EQ.HOLD) GO TO 130 120 RACUM = H/HOLD IRET1 = 1 GO TO 750 130 NQOLD = NQ TOLD = T RACUM = 1.0 IF (JSTART.GT.0) GO TO 250 GO TO 170 140 IF (JSTART.EQ.-1) GO TO 160 C ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL C DERIVATIVES ARE CALCULATED. NQ = 1 N3 = N N1 = N*10 N2 = N1 + 1 N4 = N**2 N5 = N1 + N N6 = N5 + 1 CALL DIFFUN(T,Y,SAVE(N2,1)) DO 150 I = 1,N N11 = N1 + I 150 Y(2,I) = SAVE(N11,1)*H HNEW = H K = 2 GO TO 100 C REPEAT LAST STEP BY RESTORING SAVED INFORMATION. 160 IF (NQ.EQ.NQOLD) JSTART = 1 T = TOLD NQ = NQOLD K = NQ + 1 GO TO 120 C SET THE COEFFICIENTS THAT DETERMINE THE ORDER AND THE METHOD C TYPE. CHECK FOR EXCESSIVE ORDER. THE LAST TWO STATEMENTS OF C THIS SECTION SET IWEVAL .GT.0 IF PW IS TO BE RE-EVALUATED C BECAUSE OF THE ORDER CHANGE, AND THEN REPEAT THE INTEGRATION C STEP IF IT HAS NOT YET BEEN DONE (IRET = 1) OR SKIP TO A FINAL C SCALING BEFORE EXIT IF IT HAS BEEN COMPLETED (IRET = 2). 170 IF (MF.EQ.0) GO TO 180 IF (NQ.GT.6) GO TO 190 GO TO (221,222,223,224,225,226),NQ 180 IF (NQ.GT.7) GO TO 190 GO TO (211,212,213,214,215,216,217),NQ 190 KFLAG = -2 RETURN C THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO THE MAXIMUM C ACCURACY PERMITTED BY THE MACHINE. THEY ARE IN THE ORDER USED.. C -1 C -1/2,-1/2 C -5/12,-3/4,-1/6 C -3/8,-11/12,-1/3,-1/24 C -251/720,-25/24,-35/72,-5/48,-1/120 C -95/288,-137/120,-5/8,-17/96,-1/40,-1/720 C -19087/60480,-49/40,-203/270,-49/192,-7/144,-7/1440,-1/5040 C -1 C -2/3,-1/3 C -6/11,-6/11,-1/11 C -12/25,-7/10,-1/5,-1/50 C -120/274,-225/274,-85/274,-15/274,-1/274 C -180/441,-58/63,-15/36,-25/252,-3/252,-1/1764 211 A(1) = -1.0 GO TO 230 212 A(1) = -0.500000000 A(3) = -0.500000000 GO TO 230 213 A(1) = -0.4166666666666667 A(3) = -0.750000000 A(4) = -0.1666666666666667 GO TO 230 214 A(1) = -0.375000000 A(3) = -0.9166666666666667 A(4) = -0.3333333333333333 A(5) = -0.04166666666666667 GO TO 230 215 A(1) = -0.3486111111111111 A(3) = -1.0416666666666667 A(4) = -0.4861111111111111 A(5) = -0.1041666666666667 A(6) = -0.008333333333333333 GO TO 230 216 A(1) = -0.3298611111111111 A(3) = -1.1416666666666667 A(4) = -0.625000000 A(5) = -0.1770833333333333 A(6) = -0.0250000000 A(7) = -0.001388888888888889 GO TO 230 217 A(1) = -0.3155919312169312 A(3) = -1.235000000 A(4) = -0.7518518518518519 A(5) = -0.2552083333333333 A(6) = -0.04861111111111111 A(7) = -0.004861111111111111 A(8) = -0.0001984126984126984 GO TO 230 221 A(1) = -1.000000000 GO TO 230 222 A(1) = -0.6666666666666667 A(3) = -0.3333333333333333 GO TO 230 223 A(1) = - 0.5454545454545455 A(3) = A(1) A(4) = -0.09090909090909091 GO TO 230 224 A(1) = -0.480000000 A(3) = -0.700000000 A(4) = -0.200000000 A(5) = -0.020000000 GO TO 230 225 A(1) = -0.437956204379562 A(3) = -0.8211678832116788 A(4) = -0.3102189781021898 A(5) = -0.05474452554744526 A(6) = -0.0036496350364963504 GO TO 230 226 A(1) = -0.4081632653061225 A(3) = -0.9206349206349206 A(4) = -0.4166666666666667 A(5) = -0.0992063492063492 A(6) = -0.0119047619047619 A(7) = -0.000566893424036282 230 K = NQ+1 IDOUB = K MTYP = (4 - MF)/2 ENQ2 = .5/FLOAT(NQ + 1) ENQ3 = .5/FLOAT(NQ + 2) ENQ1 = 0.5/FLOAT(NQ) PEPSH = EPS EUP = (PERTST(NQ,MTYP,2)*PEPSH)**2 E = (PERTST(NQ,MTYP,1)*PEPSH)**2 EDWN =(PERTST(NQ,MTYP,3)*PEPSH)**2 IF (EDWN.EQ.0) GO TO 780 BND = EPS*ENQ3/DBLE(FLOAT(N)) 240 IWEVAL = MF GO TO ( 250 , 680 ),IRET C THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY C MULTIPLYING THE SAVED INFORMATION BY THE PASCAL TRIANGLE C MATRIX. 250 T = T + H DO 260 J = 2,K DO 260 J1 = J,K J2 = K - J1 + J - 1 DO 260 I = 1,N 260 Y(J2,I) = Y(J2,I) + Y(J2+1,I) C UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. CONVERGENCE IS TESTED C BY REQUIRING CHANGES TO BE LESS THAN BND WHICH IS DEPENDENT ON C THE ERROR TEST CONSTANT. C THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE ARRAY C ERROR(I). IT IS EQUAL TO THE K-TH DERIVATIVE OF Y MULTIPLIED C BY H**K/(FACTORIAL(K-1)*A(K)), AND IS THEREFORE PROPORTIONAL C TO THE ACTUAL ERRORS TO THE LOWEST POWER OF H PRESENT. (H**K) DO 270 I = 1,N 270 ERROR(I) = 0.0 DO 430 L = 1,3 CALL DIFFUN (T,Y,SAVE(N2,1)) C IF THERE HAS BEEN A CHANGE OF ORDER OR THERE HAS BEEN TROUBLE C WITH CONVERGENCE, PW IS RE-EVALUATED PRIOR TO STARTING THE C CORRECTOR ITERATION IN THE CASE OF STIFF METHODS. IWEVAL IS C THEN SET TO -1 AS AN INDICATOR THAT IT HAS BEEN DONE. IF (IWEVAL.LT.1) GO TO 350 IF (MF.EQ.2) GO TO 310 CALL PEDERV(T,Y,PW,N3) R = A(1)*H DO 280 I = 1,N4 280 PW(I) = PW(I)*R 290 N11 = N3 + 1 N12 = N*N11 - N3 DO 300 I = 1,N12,N11 300 PW(I) = 1.0 + PW(I) IWEVAL = -1 CALL MATINV(PW,N,N3,J1) IF (J1.GT.0) GO TO 350 GO TO 440 310 DO 320 I = 1,N 320 SAVE(9,I) = Y(1,I) DO 340 J = 1,N R = EPS*DMAX1(EPS,DABS(SAVE(9,J))) Y(1,J) = Y(1,J) + R D = A(1)*H/R CALL DIFFUN(T,Y,SAVE(N6,1)) DO 330 I = 1,N N11 = I + (J-1)*N3 N12 = N5 + I N13 = N1 + I 330 PW(N11) = (SAVE(N12,1) - SAVE(N13,1))*D 340 Y(1,J) = SAVE(9,J) GO TO 290 350 IF (MF.NE.0) GO TO 370 DO 360 I = 1,N N11 = N1 + I 360 SAVE(9,I) = Y(2,I) - SAVE(N11,1)*H GO TO 410 370 DO 380 I = 1,N N11 = N5 + I N12 = N1 + I 380 SAVE(N11,1) = Y(2,I) - SAVE(N12,1)*H DO 400 I = 1,N D = 0.0 DO 390 J = 1,N N11 = I + (J-1)*N3 N12 = N5 + J 390 D = D + PW(N11)*SAVE(N12,1) 400 SAVE(9,I) = D 410 NT = N DO 420 I = 1,N Y(1,I) = Y(1,I) + A(1)*SAVE(9,I) Y(2,I) = Y(2,I) - SAVE(9,I) ERROR(I) = ERROR(I) + SAVE(9,I) IF (DABS(SAVE(9,I)).LE.(BND*YMAX(I))) NT = NT - 1 420 CONTINUE IF (NT.LE.0) GO TO 490 430 CONTINUE C THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. VARIOUS C POSSIBILITIES ARE CHECKED FOR. IF H IS ALREADY HMIN AND C THIS IS EITHER ADAMS METHOD OR THE STIFF METHOD IN WHICH THE C MATRIX PW HAS ALREADY BEEN RE-EVALUATED, A NO CONVERGENCE EXIT C IS TAKEN. OTHERWISE THE MATRIX PW IS RE-EVALUATED AND/OR THE C STEP IS REDUCED TO TRY AND GET CONVERGENCE. 440 T = T - H IF ((H.LE.(HMIN*1.00001)).AND.((IWEVAL - MTYP).LT.-1)) GO TO 460 IF ((MF.EQ.0).OR.(IWEVAL.NE.0)) RACUM = RACUM*0.25D0 IWEVAL = MF IRET1 = 2 GO TO 750 460 KFLAG = -3 470 DO 480 I = 1,N DO 480 J = 1,K 480 Y(J,I) = SAVE(J,I) H = HOLD NQ = NQOLD JSTART = NQ RETURN C THE CORRECTOR CONVERGED AND CONTROL IS PASSED TO STATEMENT 520 C IF THE ERROR TEST IS O.K., AND TO 540 OTHERWISE. C IF THE STEP IS O.K. IT IS ACCEPTED. IF IDOUB HAS BEEN REDUCED C TO ONE, A TEST IS MADE TO SEE IF THE STEP CAN BE INCREASED C AT THE CURRENT ORDER OR BY GOING TO ONE HIGHER OR ONE LOWER. C SUCH A CHANGE IS ONLY MADE IF THE STEP CAN BE INCREASED BY AT C LEAST 1.1. IF NO CHANGE IS POSSIBLE IDOUB IS SET TO 10 TO C PREVENT FUTHER TESTING FOR 10 STEPS. C IF A CHANGE IS POSSIBLE, IT IS MADE AND IDOUB IS SET TO C NQ + 1 TO PREVENT FURTHER TESING FOR THAT NUMBER OF STEPS. C IF THE ERROR WAS TOO LARGE, THE OPTIMUM STEP SIZE FOR THIS OR C LOWER ORDER IS COMPUTED, AND THE STEP RETRIED. IF IT SHOULD C FAIL TWICE MORE IT IS AN INDICATION THAT THE DERIVATIVES THAT C HAVE ACCUMULATED IN THE Y ARRAY HAVE ERRORS OF THE WRONG ORDER C SO THE FIRST DERIVATIVES ARE RECOMPUTED AND THE ORDER IS SET C TO 1. 490 D = 0.0 DO 500 I = 1,N 500 D = D + (ERROR(I)/YMAX(I))**2 IWEVAL = 0 IF (D.GT.E) GO TO 540 IF (K.LT.3) GO TO 520 DO 510 J = 3,K DO 510 I = 1,N 510 Y(J,I) = Y(J,I) + A(J)*ERROR(I) 520 KFLAG = +1 HNEW = H IF (IDOUB.LE.1) GO TO 550 IDOUB = IDOUB - 1 IF (IDOUB.GT.1) GO TO 700 DO 530 I = 1,N 530 SAVE(10,I) = ERROR(I) GO TO 700 540 KFLAG = KFLAG - 2 IF (H.LE.(HMIN*1.00001)) GO TO 740 T = TOLD IF (KFLAG.LE.-5) GO TO 720 550 PR2 = (D/E)**ENQ2*1.2 PR3 = 1.E+20 IF ((NQ.GE.MAXDER).OR.(KFLAG.LE.-1)) GO TO 570 D = 0.0 DO 560 I = 1,N 560 D = D + ((ERROR(I) - SAVE(10,I))/YMAX(I))**2 PR3 = (D/EUP)**ENQ3*1.4 570 PR1 = 1.E+20 IF (NQ.LE.1) GO TO 590 D = 0.0 DO 580 I = 1,N 580 D = D + (Y(K,I)/YMAX(I))**2 PR1 = (D/EDWN)**ENQ1*1.3 590 CONTINUE IF (PR2.LE.PR3) GO TO 650 IF (PR3.LT.PR1) GO TO 660 600 R = 1.0/AMAX1(PR1,1.E-4) NEWQ = NQ - 1 610 IDOUB = 10 IF ((KFLAG.EQ.1).AND.(R.LT.(1.1))) GO TO 700 IF (NEWQ.LE.NQ) GO TO 630 DO 620 I = 1,N 620 Y(NEWQ+1,I) = ERROR(I)*A(K)/DBLE(FLOAT(K)) 630 K = NEWQ + 1 IF (KFLAG.EQ.1) GO TO 670 RACUM = RACUM*R IRET1 = 3 GO TO 750 640 IF (NEWQ.EQ.NQ) GO TO 250 NQ = NEWQ GO TO 170 650 IF (PR2.GT.PR1) GO TO 600 NEWQ = NQ R = 1.0/AMAX1(PR2,1.E-4) GO TO 610 660 R = 1.0/AMAX1(PR3,1.E-4) NEWQ = NQ + 1 GO TO 610 670 IRET = 2 R = DMIN1(R,HMAX/DABS(H)) H = H*R HNEW = H IF (NQ.EQ.NEWQ) GO TO 680 NQ = NEWQ GO TO 170 680 R1 = 1.0 DO 690 J = 2,K R1 = R1*R DO 690 I = 1,N 690 Y(J,I) = Y(J,I)*R1 IDOUB = K 700 DO 710 I = 1,N 710 YMAX(I) = DMAX1(YMAX(I),DABS(Y(1,I))) JSTART = NQ RETURN 720 IF (NQ.EQ.1) GO TO 780 CALL DIFFUN (T,Y,SAVE(N2,1)) R = H/HOLD DO 730 I = 1,N Y(1,I) = SAVE(1,I) N11 = N1 + I SAVE(2,I) = HOLD*SAVE(N11,1) 730 Y(2,I) = SAVE(2,I)*R NQ = 1 KFLAG = 1 GO TO 170 740 KFLAG = -1 HNEW = H JSTART = NQ RETURN C THIS SECTION SCALES ALL VARIABLES CONNECTED WITH H AND RETURNS C TO THE ENTERING SECTION. 750 RACUM = DMAX1(DABS(HMIN/HOLD),RACUM) RACUM = DMIN1(RACUM,DABS(HMAX/HOLD)) R1 = 1.0 DO 760 J = 2,K R1 = R1*RACUM DO 760 I = 1,N 760 Y(J,I) = SAVE(J,I)*R1 H = HOLD*RACUM DO 770 I = 1,N 770 Y(1,I) = SAVE(1,I) IDOUB = K GO TO ( 130 , 250 , 640 ),IRET1 780 KFLAG = -4 GO TO 470 END