DIFSUB: PROC(N,T,Y,H,HMIN,HMAX,EPS,NQO,YMAX,ERROR,KFLAG,JSTART,MAXDER, DIFFUN) REORDER; /* */ /* THE PARAMETERS TO THE SUBROUTINE DIFSUB HAVE THE FOLLOWING */ /* MEANINGS.. */ /* */ /* N THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. N */ /* MAY BE DECREASED ON LATER CALLS IF THE NUMBER OF */ /* ACTIVE EQUATIONS REDUCES, BUT IT MUST NOT BE */ /* INCREASED WITHOUT CALLING WITH JSTART = 0. */ /* T THE INDEPENDENT VARIABLE */ /* Y AN 8 BY N ARRAY CONTAINING THE DEPENDENT VARIABLES AND */ /* THEIR SCALED DERIVATIVES. Y(J,I) CONTAINS THE J-TH */ /* DERIVATIVE OF Y(I) SCALED BY H**J/FACTORIAL(J) WHERE */ /* H IS THE CURRENT STEP SIZE. ONLY Y(0,I) NEED BE */ /* PROVIDED BY THE CALLING PROGRAM ON THE FIRST ENTRY. */ /* IF IT IS DESIRED TO INTERPOLATE TO NON-MESH POINTS */ /* THESE VALUES CAN BE USED. IF THE CURRENT STEP SIZE IS */ /* H AND THE VALUE AT T + E IS NEEDED, FORM S = E/H, AND */ /* THEN COMPUTE */ /* NQO M/ /* Y(I)(T+E) = SUM Y(J,I)*S**J */ /* J=O */ /* H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. H MAY */ /* BE ADJUSTED UP OR DOWN BY DIFSUB IN ORDER TO */ /* ACHIEVE AN ECONOMICAL INTEGRATION. HOWEVER, IF THE H */ /* PROVIDED BY THE USER DOES NOT CAUSE A LARGER ERROR */ /* THAN REQUESTED, IT WILL BE USED. TO SAVE COMPUTER */ /* TIME, THE USER IS ADVISED TO USE A FAIRLY SMALL STEP */ /* FOR THE FIRST CALL. IT WILL BE AUTOMATICALLY */ /* INCREASED LATER. */ /* HMIN THE MINIMUM STEP SIZE THAT WILL BE USED FOR THE */ /* INTEGRATION. NOTE THAT ON STARTING THIS MUST BE MUCH */ /* SMALLER THAN THE AVERAGE H EXPECTED SINCE A FIRST */ /* ORDER METHOD IS USED INITIALLY. */ /* HMAX THE MAXIMUM SIZE TO WHICH THE STEP WILL BE INCREASED */ /* EPS THE ERROR TEST CONSTANT. SINGLE STEP ERROR ESTIMATES */ /* DIVIDED BY YMAX(I) MUST BE LESS THAN THIS IN THE */ /* EUCLIDEAN NORM. THE STEP AND/OR ORDER IS ADJUSTED TO */ /* ACHIEVE THIS. */ /* NQO AN OUTPUT PARAMETER THAT IS SET TO NQ, THE CURRENT */ /* ORDER OF THE METHOD AT EXIT. NQ IS ALSO THE ORDER OF */ /* THE MAXIMUM DERIVATIVE AVAILABLE IN Y. */ /* YMAX AN ARRAY OF N ELEMENTS WHICH CONTAINS THE MAXIMUM OF */ /* EACH Y SEEN SO FAR. IT SHOULD NORMALLY BE SET TO 1 IN */ /* EACH COMPONENT BEFORE THE FIRST ENTRY. (SEE THE */ /* DESCRIPTION OF EPS.) */ /* ERROR AN ARRAY OF N ELEMENTS WHICH CONTAINS THE ESTIMATED ONE */ /* STEP ERROR IN EACH COMPONENT. */ /* KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. */ /* +1 THE STEP WAS SUCCESFUL */ /* -1 THE STEP WAS TAKEN WITH H = HMIN, BUT THE */ /* REQUESTED ERROR WAS NOT ACHIEVED. */ /* -2 THE MAXIMUM ORDER SPECIFIED WAS FOUND TO BE */ /* GREATER THAN 7. */ /* -3 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED */ /* FOR H GREATER THAN HMIN. */ /* -4 THE REQUESTED ERROR IS SMALLER THAN CAN BE */ /* HANDLED FOR THIS PROBLEM. */ /* JSTART AN INPUT INDICATOR WITH THE FOLLOWING MEANINGS.. */ /* -1 REPEAT THE LAST STEP WITH A NEW H */ /* O PERFORM THE FIRST STEP. THE FIRST STEP MUST */ /* BE DONE WITH THIS VALUE OF JSTART SO THAT */ /* THE SUBROUTINE CAN INITIALIZE ITSELF. */ /* +1 TAKE A NEW STEP CONTINUING FROM THE LAST. */ /* MAXDER THE MAXIMUM DERIVATIVE THAT SHOULD BE USED IN THE */ /* METHOD. SINCE THE ORDER IS EQUAL TO THE HIGHEST */ /* DERIVATIVE USED, THIS RESTRICTS THE ORDER. IT MUST BE */ /* LESS THAN 8 FOR ADAMS METHOD, THE METHOD USED IN */ /* DIFSUB. */ /* DIFFUN IS THE PROCEDURE SUPPLIED BY THE USER TO EVALUATE */ /* THE EQUATIONS. THE FIRST PARAMETER IS THE VALUE OF THE */ /* INDEPENDENT VARIABLE. THE SECOND PARAMETER IS Y AS */ /* DESCRIBED ABOVE. THE COMPUTED VALUES ARE RETURNED IN */ /* THE THIRD PARAMETER. */ /* */ DFT RANGE(A:H,Q:Z) BIN FLOAT VALUE (BIN FLOAT(53)); DCL (LBY,UBY) BIN FIXED(15,0); DCL Y(*,*), ERROR(*),YMAX(*),SAVE(0:7,LBY:UBY) CTL,A(8) STATIC, SAVE9(LBY:UBY) CTL, SAVE10(LBY:UBY) CTL, SAVE11(LBY:UBY) CTL, (D,E,R,R1,R2,BND,EUP,EDWN,ENQ1,ENQ2,ENQ3,HNEW,HOLD,TOLD,RACUM,K,NQ, NQOLD,IDOUB,PEPSH,NEWQ,NT,PR1,PR2,PR3,IRET,IRET1)STATIC, /* */ /* THE COEFFICIENTS IN PERTST ARE USED IN SELECTING THE STEP AND */ /* ORDER, THEREFORE ONLY ABOUT ONE PERCENT ACCURACY IS NEEDED. */ /* */ PERTST (7,3) STATIC INIT (2,12,1,12,24,1,24,37.89,2,37.89,53.33, 1,53.33,70.08,.3157,70.08,87.97,.07407,87.97,1,0.0139), L2(3) LABEL INIT(S130,S250,S640), L1(7) LABEL INIT(S211,S212,S213,S214,S215,S216,S217); DCL DIFFUN ENTRY; LBY=LBOUND(Y,2); UBY=LBY+N-1; /***************** ******************/ /* /* THE FOLLOWING CHANGES ARE BASED ON THE ARTICLE "LIMITING PRECISION /* IN DIFFERENTIAL EQUATION SOLVERS" BY L.F.SHAMPINE, MATH. OF /* COMPUTATION, VOL.28, JAN.1974, P.141. /* A CHANGE IN THE VALUE OF EPS WILL BE EFFECTIVE REGARDLESS OF THE /* VALUE OF JSTART. H AND HMIN ARE FORCED TO BE LARGE ENOUGH SO THAT /* T+H WILL DIFFER SIGNIFICANTLY FROM T. HMAX>=HMIN IS ENFORCED. /* EPS WILL BE MODIFIED IF NECESSARY TO PREVENT AN ATTEMPT TO /* COMPUTE A SOLUTION WHICH IS MORE ACCURATE THAN CAN BE OBTAINED /* WITH THE LIMITED PRECISION OF THE ARITHMETIC. */ /* */ KFLAG=1; IF H=0 THEN GO TO E1; HTEST=8E-16*T; IF HTEST^=0 THEN DO; IF ABS(HMIN)0 ! (JSTART<0 & NQ=NQOLD) THEN DO; PEPSH=EPS; EUP=(PERTST(NQ,2)*PEPSH)**2; E=(PERTST(NQ,1)*PEPSH)**2; EDWN=(PERTST(NQ,3)*PEPSH)**2; BND=EPS*ENQ3/N; END; /* */ /***************** ******************/ A(2)=-1; IRET=1; IF JSTART=0 THEN DO; IF ALLOCN(SAVE) THEN FREE SAVE,SAVE9,SAVE10,SAVE11; ALLOC SAVE,SAVE9,SAVE10,SAVE11; END; IF JSTART<=0 THEN GOTO S140; /* */ /* BEGIN BY SAVING INFORMATION FOR POSSIBLE RESTARTS AND CHANGING */ /* H BY THE FACTOR R IF THE CALLER HAS CHANGED H. ALL VARIABLES */ /* DEPENDENT ON H MUST ALSO BE CHANGED. */ /* E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS TO */ /* TEST FOR INCREASING THE ORDER, EDWN FOR DECREASING THE ORDER, */ /* HNEW IS THE STEP SIZE THAT WAS USED ON THE LAST CALL. */ /* */ S100: DO I=LBY TO UBY; DO J=0 TO K; SAVE(J,I)=Y(J,I); END; END; HOLD=HNEW; IF H=HOLD THEN GOTO S130; S120: RACUM=H/HOLD; IRET1=1; GOTO S750; S130: NQOLD=NQ; TOLD=T; RACUM=1.0; IF JSTART>0 THEN GOTO S250; GOTO S170; S140: IF JSTART=-1 THEN GOTO S160; /* */ /* ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL */ /* DERIVATIVES ARE CALCULATED. */ /* */ NQ=1; CALL DIFFUN(T,Y,SAVE11(*)); DO I=LBY TO UBY; Y(1,I)=SAVE11(I)*H; END; HNEW=H; K=1; GOTO S100; /* */ /* REPEAT LAST STEP BY RESTORING SAVED INFORMATION. */ /* */ S160: IF NQ=NQOLD THEN JSTART=1; T=TOLD; NQ=NQOLD; K=NQ; GOTO S120; /* */ /* SET THE COEFFICIENTS THAT DETERMINE THE ORDER AND THE METHOD */ /* TYPE. CHECK FOR EXCESSIVE ORDER. THE LAST STATEMENT OF THIS */ /* SECTION BRANCHES FOR A FINAL SCALING BEFORE EXIT IF INTEGRATION */ /* HAS BEEN COMPLETED. IF NOT, THE INTEGRATION IS REPEATED. */ /* */ S170: IF NQ>7 THEN DO; KFLAG=-2; NQO=0; GOTO E1; END; GOTO L1(NQ); /* */ /* THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO THE MAXIMUM */ /* ACCURACY PERMITTED BY THE MACHINE. THEY ARE IN THE ORDER USED.. */ /* */ /* -1 */ /* -1/2,-1/2 */ /* -5/12,-3/4,-1/6 */ /* -3/8,-11/12,-1/3,-1/24 */ /* -251/720,-25/24,-35/72,-5/48,-1/120 */ /* -95/288,-137/120,-5/8,-17/96,-1/40,-1/720 */ /* -19087/60480,-49/40,-203/270,-49/192,-7/144,-7/1440,-1/5040 */ /* */ S211: A(1)=-1.0; GOTO S230; S212: A(1)= -.500000000000000; A(3)= -.500000000000000; GOTO S230; S213: A(1)=-.4166666666666667E0; A(3)=-.75E0; A(4)=-.1666666666666667E0; GO TO S230; S214: A(1)=-.3750000E0; A(3)=-.9166666666666667E0; A(4)=-.3333333333333333E0; A(5)=-.4166666666666667E-1; GO TO S230; S215: A(1)=-.3486111111111111E0; A(3)=-1.041666666666667E0; A(4)=-.4861111111111111E0; A(5)=-.1041666666666667E0; A(6)=-.8333333333333333E-2; GO TO S230; S216: A(1)=-.3298611111111111E0; A(3)=-1.141666666666667E0; A(4)=-.6250000E0; A(5)=-.1770833333333333E0; A(6)=-.02500000E0; A(7)=-.1388888888888889E-2; GO TO S230; S217: A(1)=-.3155919312169312E0; A(3)=-1.2250000E0; A(4)=-.7518518518518519E0; A(5)=-.2552083333333333E0; A(6)=-.4861111111111111E-1; A(7)=-.4861111111111111E-2; A(8)=-.1984126984126984E-3; S230: K=NQ; IDOUB=K+1; ENQ2=.5000000E0/(NQ+1); ENQ3=.5000000E0/(NQ+2); ENQ1=.5000000E0/NQ; PEPSH=EPS; EUP=(PERTST(NQ,2) *PEPSH)**2; E=(PERTST(NQ,1) *PEPSH)**2; EDWN=(PERTST(NQ,3) *PEPSH)**2; IF EDWN=0.0 THEN GOTO S780; BND=EPS*ENQ3/N; S240: IF IRET=2 THEN GOTO S680; /* */ /* THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY */ /* MULTIPLYING THE SAVED INFORMATION BY THE PASCAL TRIANGLE MATRIX. */ /* */ S250: T=T+H; DO J=1 TO K; DO J1=J TO K; J2=K-J1+J-1; DO I=LBY TO UBY; Y(J2,I)=Y(J2,I)+Y(J2+1,I); END; END; END; /* */ /* UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. CONVERGENCE IS TESTED */ /* BY REQUIRING CHANGES TO BE LESS THAN BND WHICH IS DEPENDENT ON */ /* THE ERROR TEST CONSTANT. */ /* THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE ARRAY */ /* ERROR(I). IT IS EQUAL TO THE (K+1)TH DERIVATIVE OF Y MULTIPLIED */ /* BY H**(K+1)/(FACTORIAL(K)*A(K+1)), AND IS THEREFORE */ /* PROPORTIONAL TO THE ACTUAL ERRORS TO THE LOWEST POWER OF H */ /* PRESENT. (H**(K+1)) */ /* */ ERROR=0.0; DO L=1 TO 3; CALL DIFFUN(T,Y,SAVE11(*)); S350: DO I=LBY TO UBY; SAVE9(I)=Y(1,I)-SAVE11(I)*H; END; S410:NT=N; DO I=LBY TO UBY; Y(0,I)=Y(0,I)+A(1)*SAVE9(I); Y(1,I)=Y(1,I)-SAVE9(I); ERROR(I)=ERROR(I)+SAVE9(I); IF ABS(SAVE9(I))<=BND*YMAX(I) THEN NT=NT-1; S420:END; IF NT<=0 THEN GOTO S490; S430:END; /* */ /* THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. IF H IS */ /* ALREADY HMIN, A NO CONVERGENCE EXIT IS TAKEN. OTHERWISE THE */ /* STEP IS REDUCED TO TRY AND GET CONVERGENCE. */ /* */ S440:T=T-H; IF (ABS(H)<=ABS(HMIN)*1.00001) THEN GOTO S460; RACUM=RACUM*.2500; IRET1=2; GOTO S750; S460: KFLAG=-3; S470: DO I=LBY TO UBY; DO J=0 TO K; Y(J,I)=SAVE(J,I); END; END; H=HOLD; NQ=NQOLD; NQO=NQ; GOTO E1; /* */ /* THE CORRECTOR CONVERGED AND CONTROL IS PASSED TO STATEMENT S520 */ /* IF THE ERROR TEST IS O.K., AND TO S540 OTHERWISE. IF THE STEP */ /* IS O.K. IT IS ACCEPTED. IF IDOUB HAS BEEN REDUCED TO ONE, A */ /* TEST IS MADE TO SEE IF THE STEP CAN BE INCREASED AT THE CURRENT */ /* ORDER OR BY GOING TO ONE HIGHER OR ONE LOWER. SUCH A CHANGE IS */ /* ONLY MADE IF THE STEP CAN BE INCREASED BY AT LEAST 1.1. IF NO */ /* CHANGE IS POSSIBLE IDOUB IS SET TO 10 TO PREVENT FURTHER TESTING */ /* FOR 10 STEPS. */ /* IF A CHANGE IS POSSIBLE, IT IS MADE AND IDOUB IS SET TO NQ + 1 */ /* TO PREVENT FURTHER TESTING FOR THAT NUMBER OF STEPS. IF THE */ /* ERROR WAS TOO LARGE, THE OPTIMUM STEP SIZE FOR THIS OR LOWER */ /* ORDER IS COMPUTED, AND THE STEP RETRIED. IF IT SHOULD FAIL */ /* TWICE MORE IT IS AN INDICATION THAT THE DERIVATIVES THAT HAVE */ /* ACCUMULATED IN THE Y ARRAY HAVE ERRORS OF THE WRONG ORDER SO THE */ /* FIRST DERIVATIVES ARE RECOMPUTED AND THE ORDER IS SET TO 1. */ /* */ S490: D=0.0; DO I=LBY TO UBY; D=D+(ERROR(I)/YMAX(I))**2; END; /* */ /* THE STEP IS NOT ACCEPTED IF THE TRUNCATION ERROR APPROXIMATION */ /* IS GREATER THAN THE MAXIMUM ERROR ALLOWED. */ /* */ IF D>E THEN GOTO S540; DO J=2 TO K; DO I=LBY TO UBY; Y(J,I)=Y(J,I)+A(J+1)*ERROR(I); END; END; /* */ /* THE STEP IS ACCEPTED. IF IDOUB IS LESS THAN OR EQUAL TO 1, A */ /* BRANCH TO TEST THE POSSIBILITY OF INCREASING H IS TAKEN. */ /* OTHERWISE, CONTROL IS PASSED TO THE CALLING PROGRAM. */ /* */ S520: KFLAG=1; HNEW=H; IF IDOUB<=1 THEN GOTO S550; IDOUB=IDOUB-1; IF IDOUB>1 THEN GOTO S700; DO I=LBY TO UBY; SAVE10(I)=ERROR(I); END; GOTO S700; /* */ /* THE STEP IS NOT ACCEPTED. IF H IS LESS THAN H-MINIMUM, A RETURN */ /* CODE OF -1 IS ASSIGNED TO KFLAG AND CONTROL IS RETURNED TO THE */ /* CALLING PROGRAM */ /* */ S540: KFLAG=KFLAG-2; IF (ABS(H)<=ABS(HMIN)*1.00001) THEN GOTO S740; T=TOLD; IF KFLAG<=-5 THEN GOTO S720; /* */ /* PR2 IS THE STEP CHANGE FACTOR FOR DECREASING THE STEP SIZE FOR */ /* THE NEXT STEP (OR TO REPEAT IF THIS STEP FAILED). THE FACTOR OF */ /* 1.2 CAUSES THE STEP CHOSEN TO BE SMALLER THAN NEEDED TO ALLOW FOR */ /* NEGLECTED TERMS IN THE ERROR VECTOR. */ /* */ S550: PR2=(D/E)**ENQ2*1.2; PR3=1.E20; IF NQ>=MAXDER ! KFLAG<=-1 THEN GOTO S570; D=0.0; DO I=LBY TO UBY; D=D+((ERROR(I)-SAVE10(I))/YMAX(I))**2; END; /* */ /* PR3 IS THE STEP CHANGE FACTOR FOR INCREASING THE STEP SIZE. THE */ /* FACTOR 1.4 BIASES THE METHOD IN FAVOR OF CHOOSING THE SAME ORDER */ /* (NQ) OR AN ORDER ONE LESS (NQ-1). */ /* */ PR3=(D/EUP)**ENQ3*1.4; S570: PR1=1.E20; IF NQ<=1 THEN GOTO S590; D=0.0; DO I=LBY TO UBY; D=D+(Y(K,I)/YMAX(I))**2; END; /* */ /* THE ORDER MAY BE LOWERED. PR1 WILL BE USED AS AN ESTIMATE FOR */ /* THE STEP CHANGE FACTOR. THE FACTOR 1.3 BIASES THE METHOD IN */ /* FAVOR OF NOT CHANGING THE ORDER. */ /* */ PR1=(D/EDWN)**ENQ1*1.3; S590: IF PR2<=PR3 THEN GOTO S650; IF PR3PR1 THEN GOTO S600; NEWQ=NQ; R=1.0/MAX(PR2,1.E-4); GOTO S610; S660: R=1.0/MAX(PR3,1.E-4); NEWQ=NQ+1; GOTO S610; S670: IRET=2; R=MIN(R,ABS(HMAX/H)); H=H*R; HNEW=H; IF NQ=NEWQ THEN GOTO S680; NQ=NEWQ; GOTO S170; /* */ /* WHEN THE ORDER IS INCREASED FROM K TO K + 1, THE SCALED (K + 1) */ /* DERIVATIVE OF Y IS APPENDED TO VECTOR Y. */ /* */ S680: R1=1.0; DO J=1 TO K; R1=R1*R; DO I=LBY TO UBY; Y(J,I)=Y(J,I)*R1; END; END; IDOUB=K+1; /* */ /* THE MAXIMUM ACCEPTED Y(0,I) VALUE COMPUTED SO FAR IS PLACED IN */ /* YMAX(I). */ /* */ S700: DO I=LBY TO UBY; YMAX(I)=MAX(YMAX(I),ABS(Y(0,I))); END; NQO=NQ; GOTO E1; S720: IF NQ=1 THEN GOTO S780; CALL DIFFUN(T,Y,SAVE11(*)); R=H/HOLD; DO I=LBY TO UBY; Y(0,I)=SAVE(0,I); SAVE(1,I)=HOLD*SAVE11(I); Y(1,I)=SAVE(1,I)*R; END; KFLAG,NQ=1; GOTO S170; S740: KFLAG=-1; HNEW=H; NQO=NQ; GOTO E1; /* */ /* THIS SECTION SCALES ALL VARIABLES CONNECTED WITH H AND RETURNS */ /* TO THE ENTERING SECTION. */ /* */ S750: RACUM=MAX(ABS(HMIN/HOLD),RACUM); RACUM=MIN(RACUM,ABS(HMAX/HOLD)); R1=1.0; DO J=1 TO K; R1=R1*RACUM; DO I=LBY TO UBY; Y(J,I)=SAVE(J,I)*R1; END; END; H=HOLD*RACUM; DO I=LBY TO UBY; Y(0,I)=SAVE(0,I); END; IDOUB=K+1; GOTO L2(IRET1); S780: KFLAG=-4; GOTO S470; E1: END DIFSUB; STIFDIF:PROC(N,T,Y,H,HMIN,HMAX,EPS,MF,YMAX,ERROR,KFLAG,JSTART,MAXDER, NQO,DIFFUN,PEDERV) REORDER; /* */ /* THE PARAMETERS TO PROCEDURE STIFDIF HAVE THE FOLLOWING MEANINGS: */ /* */ /* N THE NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS. N */ /* MAY BE DECREASED ON LATER CALLS IF THE NUMBER OF */ /* ACTIVE EQUATIONS REDUCES, BUT IT MUST NOT BE */ /* INCREASED WITHOUT CALLING WITH JSTART = 0. */ /* T THE INDEPENDENT VARIABLE */ /* Y AN 8 BY N ARRAY CONTAINING THE DEPENDENT VARIABLES AND */ /* THEIR SCALED DERIVATIVES. Y(J,I) CONTAINS THE J-TH */ /* DERIVATIVE OF Y(I) SCALED BY H**J/FACTORIAL(J) WHERE */ /* H IS THE CURRENT STEP SIZE. ONLY Y(0,I) NEED BE */ /* PROVIDED BY THE CALLING PROGRAM ON THE FIRST ENTRY. */ /* IF IT IS DESIRED TO INTERPOLATE TO NON-MESH POINTS */ /* THESE VALUES CAN BE USED. IF THE CURRENT STEP SIZE IS */ /* H AND THE VALUE AT T + E IS NEEDED, FORM S = E/H, AND */ /* THEN COMPUTE */ /* NQO M/ /* Y(I)(T+E) = SUM Y(J,I)*S**J */ /* J=O */ /* H THE STEP SIZE TO BE ATTEMPTED ON THE NEXT STEP. H MAY */ /* BE ADJUSTED UP OR DOWN BY DIFSUB IN ORDER TO */ /* ACHIEVE AN ECONOMICAL INTEGRATION. HOWEVER, IF THE H */ /* PROVIDED BY THE USER DOES NOT CAUSE A LARGER ERROR */ /* THAN REQUESTED, IT WILL BE USED. TO SAVE COMPUTER */ /* TIME, THE USER IS ADVISED TO USE A FAIRLY SMALL STEP */ /* FOR THE FIRST CALL. IT WILL BE AUTOMATICALLY */ /* INCREASED LATER. */ /* HMIN THE MINIMUM STEP SIZE THAT WILL BE USED FOR THE */ /* INTEGRATION. NOTE THAT ON STARTING THIS MUST BE MUCH */ /* SMALLER THAN THE AVERAGE H EXPECTED SINCE A FIRST */ /* ORDER METHOD IS USED INITIALLY. */ /* HMAX THE MAXIMUM SIZE TO WHICH THE STEP WILL BE INCREASED */ /* EPS THE ERROR TEST CONSTANT. SINGLE STEP ERROR ESTIMATES */ /* DIVIDED BY YMAX(I) MUST BE LESS THAN THIS IN THE */ /* EUCLIDEAN NORM. THE STEP AND/OR ORDER IS ADJUSTED TO */ /* ACHIEVE THIS. */ /* MF THE METHOD INDICATOR. THE FOLLOWING ARE ALLOWED.. */ /* 1 A MULTI-STEP METHOD SUITABLE FOR STIFF SYSTEMS IS USED. */ /* IT WILL ALSO WORK FOR NON STIFF SYSTEMS. HOWEVER THE */ /* USER MUST PROVIDE A PROCEDURE PEDERV WHICH EVALUATES */ /* THE PARTIAL DERIVATIVES OF THE DIFFERENTIAL EQUATIONS */ /* WITH RESPECT TO THE Y'S. */ /* 2 THE SAME AS CASE 1, EXCEPT THAT THIS PROCEDURE COMPUTES */ /* THE PARTIAL DERIVATIVES BY NUMERICAL DIFFERENCING */ /* OF THE DERIVATIVES. HENCE PEDERV IS NOT CALLED. */ /* YMAX AN ARRAY OF N ELEMENTS WHICH CONTAINS THE MAXIMUM OF */ /* EACH Y SEEN SO FAR. IT SHOULD NORMALLY BE SET TO 1 IN */ /* EACH COMPONENT BEFORE THE FIRST ENTRY. (SEE THE */ /* DESCRIPTION OF EPS.) */ /* ERROR AN ARRAY OF N ELEMENTS WHICH CONTAINS THE ESTIMATED ONE */ /* STEP ERROR IN EACH COMPONENT. */ /* KFLAG A COMPLETION CODE WITH THE FOLLOWING MEANINGS.. */ /* +1 THE STEP WAS SUCCESFUL */ /* -1 THE STEP WAS TAKEN WITH H = HMIN, BUT THE */ /* REQUESTED ERROR WAS NOT ACHIEVED. */ /* -2 THE MAXIMUM ORDER SPECIFIED WAS FOUND TO BE */ /* GREATER THAN 6. */ /* -3 CORRECTOR CONVERGENCE COULD NOT BE ACHIEVED */ /* FOR H GREATER THAN HMIN. */ /* -4 THE REQUESTED ERROR IS SMALLER THAN CAN BE */ /* HANDLED FOR THIS PROBLEM. */ /* JSTART AN INPUT INDICATOR WITH THE FOLLOWING MEANINGS.. */ /* -1 REPEAT THE LAST STEP WITH A NEW H */ /* O PERFORM THE FIRST STEP. THE FIRST STEP MUST */ /* BE DONE WITH THIS VALUE OF JSTART SO THAT */ /* THE SUBROUTINE CAN INITIALIZE ITSELF. */ /* +1 TAKE A NEW STEP CONTINUING FROM THE LAST. */ /* MAXDER THE MAXIMUM DERIVATIVE THAT SHOULD BE USED IN THE */ /* METHOD. SINCE THE ORDER IS EQUAL TO THE HIGHEST */ /* DERIVATIVE USED, THIS RESTRICTS THE ORDER. IT MUST BE */ /* LESS THAN 7 FOR THE ADAM'S METHOD USED IN STIFDIF. */ /* NQO AN OUTPUT PARAMETER THAT IS SET TO NQ, THE CURRENT */ /* ORDER OF THE METHOD AT EXIT. NQ IS ALSO THE ORDER OF */ /* THE MAXIMUM DERIVATIVE AVAILABLE IN Y. */ /* DIFFUN IS THE PROCEDURE SUPPLIED BY THE USER TO EVALUATE */ /* THE EQUATIONS. THE FIRST PARAMETER IS THE VALUE OF THE */ /* INDEPENDENT VARIABLE. THE SECOND PARAMETER IS Y AS */ /* DESCRIBED ABOVE. THE COMPUTED VALUES ARE RETURNED IN */ /* THE THIRD PARAMETER. */ /* PEDERV THE NAME OF THE PROCEDURE SUPPLIED BY THE USER TO */ /* EVALUATE THE PARTIAL DERIVATIVES OF THE DIFFERENTIAL */ /* EQUATIONS. */ /* */ DFT RANGE(A:H,Q:Z) BIN FLOAT VALUE (BIN FLOAT(53)); DCL Y(*,*), ERROR(*),YMAX(*),A(8) STATIC,UBY BIN FIXED, SAVE(0:7,LBY:UBY) CTL,SAVE10(LBY:UBY) CTL, SAVE11(LBY:UBY) CTL,SAVE12(LBY:UBY) CTL,PW(LBY:UBY,LBY:UBY) CTL, PW1(LBY:UBY,LBY:UBY) CTL, INT(LBY:UBY) BIN FIXED (31) CTL, /* */ /* THE COEFFICIENTS IN PERTST ARE USED IN SELECTING THE STEP AND */ /* ORDER, THEREFORE ONLY ABOUT ONE PERCENT ACCURACY IS NEEDED. */ /* */ PERTST (6,3) STATIC INIT (2,3,1,4.5,6,1,7.333,9.167,.5,10.42, 12.5,.1667,13.7,15.98,.04133,17.15,1,.008267), (D,E,R,R1,BND,EUP,EDWN,ENQ1,ENQ2,ENQ3,HNEW,HOLD,TOLD,RACUM,K,NQ, NQOLD,IDOUB,PEPSH,IWEVAL,NEWQ,NT,PR1,PR2,PR3,IRET,IRET1) STATIC, (LL,D2) BIN FIXED(31), L2(3) LABEL INIT(S130,S250,S640), L1(6) LABEL INIT(S221,S222,S223,S224,S225,S226); DCL UNSLVE EXT ENTRY (BIN FIXED (31),BIN FIXED (31),,,,,,,,,), UNSDET EXT ENTRY (BIN FIXED (31),,,,,,); DCL (B,BB,X) (LBY:UBY,1) BIN FLOAT (21) CTL; DCL (DIFFUN,PEDERV) ENTRY; LBY=LBOUND(Y,2); UBY=LBY+N-1; /***************** ******************/ /* /* THE FOLLOWING CHANGES ARE BASED ON THE ARTICLE "LIMITING PRECISION /* IN DIFFERENTIAL EQUATION SOLVERS" BY L.F.SHAMPINE, MATH. OF /* COMPUTATION, VOL.28, JAN.1974, P.141. /* A CHANGE IN THE VALUE OF EPS WILL BE EFFECTIVE REGARDLESS OF THE /* VALUE OF JSTART. H AND HMIN ARE FORCED TO BE LARGE ENOUGH SO THAT /* T+H WILL DIFFER SIGNIFICANTLY FROM T. HMAX>=HMIN IS ENFORCED. /* EPS WILL BE MODIFIED IF NECESSARY TO PREVENT AN ATTEMPT TO /* COMPUTE A SOLUTION WHICH IS MORE ACCURATE THAN CAN BE OBTAINED /* WITH THE LIMITED PRECISION OF THE ARITHMETIC. */ /* */ KFLAG=1; IF H=0 THEN GO TO E1; HTEST=8E-16*T; IF HTEST^=0 THEN DO; IF ABS(HMIN)0 ! (JSTART<0 & NQ=NQOLD) THEN DO; PEPSH=EPS; EUP=(PERTST(NQ,2)*PEPSH)**2; E=(PERTST(NQ,1)*PEPSH)**2; EDWN=(PERTST(NQ,3)*PEPSH)**2; BND=EPS*ENQ3/N; END; /* */ /***************** ******************/ A(2)=-1; IRET=1; IF JSTART = 0 THEN DO; IF ALLOCN(SAVE) THEN FREE SAVE,SAVE10, PW,PW1,B,X,BB,SAVE11,SAVE12,INT; ALLOC SAVE,SAVE10,SAVE11,SAVE12,INT,PW,PW1,B,X,BB; END; IF JSTART<=0 THEN GOTO S140; /* */ /* BEGIN BY SAVING INFORMATION FOR POSSIBLE RESTARTS AND CHANGING */ /* H BY THE FACTOR R IF THE CALLER HAS CHANGED H. ALL VARIABLES */ /* DEPENDENT ON H MUST ALSO BE CHANGED. */ /* E IS A COMPARISON FOR ERRORS OF THE CURRENT ORDER NQ. EUP IS TO */ /* TEST FOR INCREASING THE ORDER, EDWN FOR DECREASING THE ORDER, */ /* HNEW IS THE STEP SIZE THAT WAS USED ON THE LAST CALL. */ /* */ S100: DO I=LBY TO UBY; DO J=0 TO K; SAVE(J,I)=Y(J,I); END; END; HOLD=HNEW; IF H=HOLD THEN GOTO S130; S120: RACUM=H/HOLD; IRET1=1; GOTO S750; S130: NQOLD=NQ; TOLD=T; RACUM=1.0; IF JSTART>0 THEN GOTO S250; GOTO S170; S140: IF JSTART=-1 THEN GOTO S160; /* */ /* ON THE FIRST CALL, THE ORDER IS SET TO 1 AND THE INITIAL */ /* DERIVATIVES ARE CALCULATED. */ /* */ NQ=1; CALL DIFFUN(T,Y,SAVE11(*)); DO I=LBY TO UBY; Y(1,I)=SAVE11(I)*H; END; HNEW=H; K=1; GOTO S100; /* */ /* REPEAT LAST STEP BY RESTORING SAVED INFORMATION. */ /* */ S160: IF NQ=NQOLD THEN JSTART=1; T=TOLD; NQ=NQOLD; K=NQ; GOTO S120; /* */ /* SET THE COEFFICIENTS THAT DETERMINE THE ORDER AND THE METHOD */ /* TYPE. CHECK FOR EXCESSIVE ORDER. THE LAST STATEMENT OF THIS */ /* SECTION BRANCHES FOR A FINAL SCALING BEFORE EXIT IF INTEGRATION */ /* HAS BEEN COMPLETED. IF NOT, THE INTEGRATION IS REPEATED. */ /* */ S170: IF NQ>6 THEN DO; KFLAG=-2; NQO=0; GOTO E1; END; GOTO L1(NQ); /* */ /* THE FOLLOWING COEFFICIENTS SHOULD BE DEFINED TO THE MAXIMUM */ /* ACCURACY PERMITTED BY THE MACHINE. THEY ARE IN THE ORDER USED.. */ /* */ /* -1 */ /* -2/3,-1/3 */ /* -6/11,-6/11,-1/11 */ /* -12/25,-7/10,-1/5,-1/50 */ /* -120/274,-225/274,-85/274,-15/274,-1/274 */ /* -180/441,-58/63,-15/36,-25/252,-3/252,-1/1764 */ /* */ S221: A(1)=-1.0; GOTO S230; S222: A(1)=-.6666666666666667E0; A(3)=-.3333333333333333E0; GO TO S230; S223: A(1)=-.5454545454545455E0; A(3)=A(1); A(4)=-.9090909090909091E-1; GO TO S230; S224: A(1)=-.4800000E0; A(3)=-.7000000E0; A(4)=-.20000000E0; A(5)=-.2000000E-1; GO TO S230; S225: A(1)=-.4379562043795620E0; A(3)=-.8211678832116788E0; A(4)=-.3102189781021898E0; A(5)=-.5474452554744526E-1; A(6)=-.364963503649635E-2; GO TO S230; S226: A(1)=-.4081632653061225E0; A(3)=-.9206349206349206E0; A(4)=-.4166666666666667E0; A(5)=-.9920634920634921E-1; A(6)=-.119047619047619E-1; A(7)=-.566893424036282E-3; S230: K=NQ; IDOUB=K+1; ENQ2=.5000000E0/(NQ+1); ENQ3=.5000000E0/(NQ+2); ENQ1=.5000000E0/NQ; PEPSH=EPS; EUP=(PERTST(NQ,2)*PEPSH)**2; E=(PERTST(NQ,1)*PEPSH)**2; EDWN=(PERTST(NQ,3)*PEPSH)**2; IF EDWN=0.0 THEN GOTO S780; BND=EPS*ENQ3/N; S240: IWEVAL=MF; IF IRET=2 THEN GOTO S680; /* */ /* THIS SECTION COMPUTES THE PREDICTED VALUES BY EFFECTIVELY */ /* MULTIPLYING THE SAVED INFORMATION BY THE PASCAL TRIANGLE MATRIX. */ /* */ S250: T=T+H; DO J=1 TO K; DO J1=J TO K; J2=K-J1+J-1; DO I=LBY TO UBY; Y(J2,I)=Y(J2,I)+Y(J2+1,I); END; END; END; /* */ /* UP TO 3 CORRECTOR ITERATIONS ARE TAKEN. CONVERGENCE IS TESTED */ /* BY REQUIRING CHANGES TO BE LESS THAN BND WHICH IS DEPENDENT ON */ /* THE ERROR TEST CONSTANT. */ /* THE SUM OF THE CORRECTIONS IS ACCUMULATED IN THE ARRAY */ /* ERROR(I). IT IS EQUAL TO THE (K+1)TH DERIVATIVE OF Y MULTIPLIED */ /* BY H**(K+1)/(FACTORIAL(K)*A(K+1)), AND IS THEREFORE */ /* PROPORTIONAL TO THE ACTUAL ERRORS TO THE LOWEST POWER OF H */ /* PRESENT. (H**(K+1)) */ /* */ ERROR=0.0; DO L=1 TO 3; CALL DIFFUN(T,Y,SAVE11(*)); IF IWEVAL < 1 THEN GO TO S370; IF MF=2 THEN DO; S310: DO J=LBY TO UBY; SAVE9=Y(0,J); R=MAX(ABS(SAVE9)*MAX(8E-16,EPS),EPS**2); Y(0,J)=Y(0,J)+R; D=A(1)*H/R; CALL DIFFUN(T,Y,SAVE12(*)); DO I=LBY TO UBY; PW(I,J)=(SAVE12(I)-SAVE11(I))*D; END; Y(0,J)=SAVE9; END; GOTO S290; END; CALL PEDERV(T,Y,PW); R=A(1)*H; PW=PW*R; S290: DO I=LBY TO UBY; PW(I,I)=PW(I,I)+1.0; END; PW1=PW; IWEVAL=-1; CALL UNSDET(N,.953674E-6,PW1,D1,D2,INT,S440); S370: DO I=LBY TO UBY; B(I,1)=Y(1,I)-SAVE11(I)*H; END; CALL UNSLVE (N,1,PW,PW1,INT,B,.953674E-6,X,BB,LL,S440); NT=N; DO I=LBY TO UBY; Y(0,I)=Y(0,I)+A(1)*X(I,1); Y(1,I)=Y(1,I)-X(I,1); ERROR(I)=ERROR(I)+X(I,1); IF ABS(X(I,1))<=BND*YMAX(I) THEN NT=NT-1; END; IF NT<=0 THEN GOTO S490; END; /* */ /* THE CORRECTOR ITERATION FAILED TO CONVERGE IN 3 TRIES. IF H IS */ /* ALREADY HMIN, A NO CONVERGENCE EXIT IS TAKEN. OTHERWISE THE */ /* STEP IS REDUCED TO TRY AND GET CONVERGENCE. */ /* */ S440: T=T-H; IF ABS(H)<=ABS(HMIN)*1.00001 & IWEVAL<0 THEN GO TO S460; IF IWEVAL^=0 THEN RACUM=RACUM*.2500; IWEVAL=MF; IRET1=2; GOTO S750; S460: KFLAG=-3; S470: DO I=LBY TO UBY; DO J=0 TO K; Y(J,I)=SAVE(J,I); END; END; H=HOLD; NQ=NQOLD; NQO=NQ; GOTO E1; /* */ /* THE CORRECTOR CONVERGED AND CONTROL IS PASSED TO STATEMENT S520 */ /* IF THE ERROR TEST IS O.K., AND TO S540 OTHERWISE. IF THE STEP */ /* IS O.K. IT IS ACCEPTED. IF IDOUB HAS BEEN REDUCED TO ONE, A */ /* TEST IS MADE TO SEE IF THE STEP CAN BE INCREASED AT THE CURRENT */ /* ORDER OR BY GOING TO ONE HIGHER OR ONE LOWER. SUCH A CHANGE IS */ /* ONLY MADE IF THE STEP CAN BE INCREASED BY AT LEAST 1.1. IF NO */ /* CHANGE IS POSSIBLE IDOUB IS SET TO 10 TO PREVENT FURTHER TESTING */ /* FOR 10 STEPS. */ /* IF A CHANGE IS POSSIBLE, IT IS MADE AND IDOUB IS SET TO NQ + 1 */ /* TO PREVENT FURTHER TESTING FOR THAT NUMBER OF STEPS. IF THE */ /* ERROR WAS TOO LARGE, THE OPTIMUM STEP SIZE FOR THIS OR LOWER */ /* ORDER IS COMPUTED, AND THE STEP RETRIED. IF IT SHOULD FAIL */ /* TWICE MORE IT IS AN INDICATION THAT THE DERIVATIVES THAT HAVE */ /* ACCUMULATED IN THE Y ARRAY HAVE ERRORS OF THE WRONG ORDER SO THE */ /* FIRST DERIVATIVES ARE RECOMPUTED AND THE ORDER IS SET TO 1. */ /* */ S490: D=0.0; DO I=LBY TO UBY; D=D+(ERROR(I)/YMAX(I))**2; END; /* */ /* THE STEP IS NOT ACCEPTED IF THE TRUNCATION ERROR APPROXIMATION */ /* IS GREATER THAN THE MAXIMUM ERROR ALLOWED. */ /* */ IWEVAL=0; IF D>E THEN GOTO S540; DO J=2 TO K; DO I=LBY TO UBY; Y(J,I)=Y(J,I)+A(J+1)*ERROR(I); END; END; /* */ /* THE STEP IS ACCEPTED. IF IDOUB IS LESS THAN OR EQUAL TO 1, A */ /* BRANCH TO TEST THE POSSIBILITY OF INCREASING H IS TAKEN. */ /* OTHERWISE, CONTROL IS PASSED TO THE CALLING PROGRAM. */ /* */ S520: KFLAG=1; HNEW=H; IF IDOUB<=1 THEN GOTO S550; IDOUB=IDOUB-1; IF IDOUB>1 THEN GOTO S700; DO I=LBY TO UBY; SAVE10(I)=ERROR(I); END; GOTO S700; /* */ /* THE STEP IS NOT ACCEPTED. IF H IS LESS THAN H-MINIMUM, A RETURN */ /* CODE OF -1 IS ASSIGNED TO KFLAG AND CONTROL IS RETURNED TO THE */ /* CALLING PROGRAM */ /* */ S540: KFLAG=KFLAG-2; IF ABS(H)<=ABS(HMIN)*1.00001 THEN GOTO S740; T=TOLD; IF KFLAG<=-5 THEN GOTO S720; /* */ /* PR2 IS THE STEP CHANGE FACTOR FOR DECREASING THE STEP SIZE FOR */ /* THE NEXT STEP (OR TO REPEAT IF THIS STEP FAILED). THE FACTOR OF */ /* 1.2 CAUSES THE STEP CHOSEN TO BE SMALLER THAN NEEDED TO ALLOW FOR */ /* NEGLECTED TERMS IN THE ERROR VECTOR. */ /* */ S550: PR2=(D/E)**ENQ2*1.2; PR3=1.E20; IF NQ>=MAXDER ! KFLAG<=-1 THEN GOTO S570; D=0.0; DO I=LBY TO UBY; D=D+((ERROR(I)-SAVE10(I))/YMAX(I))**2; END; /* */ /* PR3 IS THE STEP CHANGE FACTOR FOR INCREASING THE STEP SIZE. THE */ /* FACTOR 1.4 BIASES THE METHOD IN FAVOR OF CHOOSING THE SAME ORDER */ /* (NQ) OR AN ORDER ONE LESS (NQ-1). */ /* */ PR3=(D/EUP)**ENQ3*1.4; S570: PR1=1.E20; IF NQ<=1 THEN GOTO S590; D=0.0; DO I=LBY TO UBY; D=D+(Y(K,I)/YMAX(I))**2; END; /* */ /* THE ORDER MAY BE LOWERED. PR1 WILL BE USED AS AN ESTIMATE FOR */ /* THE STEP CHANGE FACTOR. THE FACTOR 1.3 BIASES THE METHOD IN */ /* FAVOR OF NOT CHANGING THE ORDER. */ /* */ PR1=(D/EDWN)**ENQ1*1.3; S590: IF PR2<=PR3 THEN GOTO S650; IF PR3PR1 THEN GOTO S600; NEWQ=NQ; R=1.0/MAX(PR2,1.E-4); GOTO S610; S660:R=1.0/MAX(PR3,1.E-4); NEWQ=NQ+1; GOTO S610; S670: IRET=2; R=MIN(R,ABS(HMAX/H)); H=H*R; HNEW=H; IF NQ=NEWQ THEN GOTO S680; NQ=NEWQ; GOTO S170; /* */ /* WHEN THE ORDER IS INCREASED FROM K TO K + 1, THE SCALED (K + 1) */ /* DERIVATIVE OF Y IS APPENDED TO VECTOR Y. */ /* */ S680: R1=1.0; DO J=1 TO K; R1=R1*R; DO I=LBY TO UBY; Y(J,I)=Y(J,I)*R1; END; END; IDOUB=K+1; /* */ /* THE MAXIMUM ACCEPTED Y(0,I) VALUE COMPUTED SO FAR IS PLACED IN */ /* YMAX(I). */ /* */ S700: DO I=LBY TO UBY; YMAX(I)=MAX(YMAX(I),ABS(Y(0,I))); END; NQO=NQ; GOTO E1; S720: IF NQ=1 THEN GOTO S780; CALL DIFFUN(T,Y,SAVE11(*)); R=H/HOLD; DO I=LBY TO UBY; Y(0,I)=SAVE(0,I); SAVE(1,I)=HOLD*SAVE11(I); Y(1,I)=SAVE(1,I)*R; END; KFLAG,NQ=1; GOTO S170; S740: KFLAG=-1; HNEW=H; NQO=NQ; GOTO E1; /* */ /* THIS SECTION SCALES ALL VARIABLES CONNECTED WITH H AND RETURNS */ /* TO THE ENTERING SECTION. */ /* */ S750: RACUM=MAX(ABS(HMIN/HOLD),RACUM); RACUM=MIN(RACUM,ABS(HMAX/HOLD)); R1=1.0; DO J=1 TO K; R1=R1*RACUM; DO I=LBY TO UBY; Y(J,I)=SAVE(J,I)*R1; END; END; H=HOLD*RACUM; DO I=LBY TO UBY; Y(0,I)=SAVE(0,I); END; IDOUB=K+1; GOTO L2(IRET1); S780: KFLAG=-4; GOTO S470; E1: END STIFDIF;