C C .................................................................. C C SUBROUTINE DHPCG C C PURPOSE C TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY GENERAL C DIFFERENTIAL EQUATIONS WITH GIVEN INITIAL VALUES. C>>>>>> THE ORIGINAL I.B.M. PROGRAM WAS MODIFIED IN TWO WAYS: C>>>>>> 1) THE INTEGRATION STEP H IS ALLOWED TO GROW LARGER C>>>>>> THAN THE INITIAL VALUE PRMT(3) . C>>>>>> 2) A PARITY CHECK WAS ADDED TO GUARANTEE THAT (X-PRMT(1))/H C>>>>>> ALWAYS REMAINS INTEGER. THUS, IF THE INTEGRATION INTERVAL C>>>>>> PRMT(2)-PRMT(1) = PRMT(3)*2**N , THEN THE INTEGRATION C>>>>>> WILL ALWAYS STOP EXACTLY AT PRMT(2) , REGARDLESS OF THE C>>>>>> VALUE OF H . (PROGRAMMER: A. YAHIL, VERSION: 4-3-79) C C USAGE C CALL DHPCG (PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) C PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT. C C DESCRIPTION OF PARAMETERS C PRMT - DOUBLE PRECISION INPUT AND OUTPUT VECTOR WITH C DIMENSION GREATER THAN OR EQUAL TO 5, WHICH C SPECIFIES THE PARAMETERS OF THE INTERVAL AND OF C ACCURACY AND WHICH SERVES FOR COMMUNICATION BETWEEN C OUTPUT SUBROUTINE (FURNISHED BY THE USER) AND C SUBROUTINE DHPCG. EXCEPT PRMT(5) THE COMPONENTS C ARE NOT DESTROYED BY SUBROUTINE DHPCG AND THEY ARE C PRMT(1)- LOWER BOUND OF THE INTERVAL (INPUT), C PRMT(2)- UPPER BOUND OF THE INTERVAL (INPUT), C PRMT(3)- INITIAL INCREMENT OF THE INDEPENDENT VARIABLE C (INPUT), C PRMT(4)- UPPER ERROR BOUND (INPUT). IF ABSOLUTE ERROR IS C GREATER THAN PRMT(4), INCREMENT GETS HALVED. C IF INCREMENT IS LESS THAN PRMT(3) AND ABSOLUTE C ERROR LESS THAN PRMT(4)/50, INCREMENT GETS DOUBLED. C THE USER MAY CHANGE PRMT(4) BY MEANS OF HIS C OUTPUT SUBROUTINE. C PRMT(5)- NO INPUT PARAMETER. SUBROUTINE DHPCG INITIALIZES C PRMT(5)=0. IF THE USER WANTS TO TERMINATE C SUBROUTINE DHPCG AT ANY OUTPUT POINT, HE HAS TO C CHANGE PRMT(5) TO NON-ZERO BY MEANS OF SUBROUTINE C OUTP. FURTHER COMPONENTS OF VECTOR PRMT ARE C FEASIBLE IF ITS DIMENSION IS DEFINED GREATER C THAN 5. HOWEVER SUBROUTINE DHPCG DOES NOT REQUIRE C AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL C FOR HANDING RESULT VALUES TO THE MAIN PROGRAM C (CALLING DHPCG) WHICH ARE OBTAINED BY SPECIAL C MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP. C Y - DOUBLE PRECISION INPUT VECTOR OF INITIAL VALUES C (DESTROYED). LATERON Y IS THE RESULTING VECTOR OF C DEPENDENT VARIABLES COMPUTED AT INTERMEDIATE C POINTS X. C DERY - DOUBLE PRECISION INPUT VECTOR OF ERROR WEIGHTS C (DESTROYED). THE SUM OF ITS COMPONENTS MUST BE C EQUAL TO 1. LATERON DERY IS THE VECTOR OF C DERIVATIVES, WHICH BELONG TO FUNCTION VALUES Y AT C INTERMEDIATE POINTS X. C NDIM - AN INPUT VALUE, WHICH SPECIFIES THE NUMBER OF C EQUATIONS IN THE SYSTEM. C IHLF - AN OUTPUT VALUE, WHICH SPECIFIES THE NUMBER OF C BISECTIONS OF THE INITIAL INCREMENT. IF IHLF GETS C GREATER THAN 10, SUBROUTINE DHPCG RETURNS WITH C ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM. C ERROR MESSAGE IHLF=12 OR IHLF=13 APPEARS IN CASE C PRMT(3)=0 OR IN CASE SIGN(PRMT(3)).NE.SIGN(PRMT(2)- C PRMT(1)) RESPECTIVELY. C FCT - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT C COMPUTES THE RIGHT HAND SIDES DERY OF THE SYSTEM C TO GIVEN VALUES OF X AND Y. ITS PARAMETER LIST C MUST BE X,Y,DERY. THE SUBROUTINE SHOULD NOT C DESTROY X AND Y. C OUTP - THE NAME OF AN EXTERNAL OUTPUT SUBROUTINE USED. C ITS PARAMETER LIST MUST BE X,Y,DERY,IHLF,NDIM,PRMT. C NONE OF THESE PARAMETERS (EXCEPT, IF NECESSARY, C PRMT(4),PRMT(5),...) SHOULD BE CHANGED BY C SUBROUTINE OUTP. IF PRMT(5) IS CHANGED TO NON-ZERO, C SUBROUTINE DHPCG IS TERMINATED. C AUX - DOUBLE PRECISION AUXILIARY STORAGE ARRAY WITH 16 C ROWS AND NDIM COLUMNS. C C REMARKS C THE PROCEDURE TERMINATES AND RETURNS TO CALLING PROGRAM, IF C (1) MORE THAN 10 BISECTIONS OF THE INITIAL INCREMENT ARE C NECESSARY TO GET SATISFACTORY ACCURACY (ERROR MESSAGE C IHLF=11), C (2) INITIAL INCREMENT IS EQUAL TO 0 OR HAS WRONG SIGN C (ERROR MESSAGES IHLF=12 OR IHLF=13), C (3) THE WHOLE INTEGRATION INTERVAL IS WORKED THROUGH, C (4) SUBROUTINE OUTP HAS CHANGED PRMT(5) TO NON-ZERO. C C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED C THE EXTERNAL SUBROUTINES FCT(X,Y,DERY) AND C OUTP(X,Y,DERY,IHLF,NDIM,PRMT) MUST BE FURNISHED BY THE USER. C C METHOD C EVALUATION IS DONE BY MEANS OF HAMMINGS MODIFIED PREDICTOR- C CORRECTOR METHOD. IT IS A FOURTH ORDER METHOD, USING 4 C PRECEEDING POINTS FOR COMPUTATION OF A NEW VECTOR Y OF THE C DEPENDENT VARIABLES. C FOURTH ORDER RUNGE-KUTTA METHOD SUGGESTED BY RALSTON IS C USED FOR ADJUSTMENT OF THE INITIAL INCREMENT AND FOR C COMPUTATION OF STARTING VALUES. C SUBROUTINE DHPCG AUTOMATICALLY ADJUSTS THE INCREMENT DURING C THE WHOLE COMPUTATION BY HALVING OR DOUBLING. C TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE C MUST BE CODED BY THE USER. C FOR REFERENCE, SEE C (1) RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL C COMPUTERS, WILEY, NEW YORK/LONDON, 1960, PP.95-109. C (2) RALSTON, RUNGE-KUTTA METHODS WITH MINIMUM ERROR BOUNDS, C MTAC, VOL.16, ISS.80 (1962), PP.431-437. C C .................................................................. C SUBROUTINE DHPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) C C implicit real*8 (a-h,o-z) real*8 :: PRMT(:),Y(:),DERY(:),AUX(:,:) real*8 ::X,H,Z,DELT c interface c c c subroutine outp(x,y,dery,ihlf,ndim,prmt) c ---------------------------------------- implicit none real*8 :: x,y(:),dery(:),prmt(:) integer :: ihlf,ndim end subroutine outp c c c c subroutine fct(x,y,dery) c ----------------------- implicit none real*8 :: x,y(:),dery(:) end subroutine fct end interface C C ---------- IPAR DECLARED REAL SO VAX CAN HANDLE SMALL STEPS C REAL*4 IPAR N=1 IHLF=0 X=PRMT(1) H=PRMT(3) PRMT(5)=0.D0 DO 1 I=1,NDIM AUX(16,I)=0.D0 AUX(15,I)=DERY(I) 1 AUX(1,I)=Y(I) IF(H*(PRMT(2)-X))3,2,4 C C ERROR RETURNS 2 IHLF=12 GOTO 4 3 IHLF=13 C C COMPUTATION OF DERY FOR STARTING VALUES 4 CALL FCT(X,Y,DERY) C C RECORDING OF STARTING VALUES CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))6,5,6 5 IF(IHLF)7,7,6 6 RETURN 7 DO 8 I=1,NDIM 8 AUX(8,I)=DERY(I) C C COMPUTATION OF AUX(2,I) ISW=1 GOTO 100 C 9 X=X+H DO 10 I=1,NDIM 10 AUX(2,I)=Y(I) C C INCREMENT H IS TESTED BY MEANS OF BISECTION 11 IHLF=IHLF+1 X=X-H DO 12 I=1,NDIM 12 AUX(4,I)=AUX(2,I) H=.5D0*H N=1 ISW=2 GOTO 100 C 13 X=X+H CALL FCT(X,Y,DERY) N=2 DO 14 I=1,NDIM AUX(2,I)=Y(I) 14 AUX(9,I)=DERY(I) ISW=3 GOTO 100 C C COMPUTATION OF TEST VALUE DELT 15 DELT=0.D0 DO 16 I=1,NDIM 16 DELT=DELT+AUX(15,I)*DABS(Y(I)-AUX(4,I)) DELT=.066666666666666667D0*DELT IF(DELT-PRMT(4))19,19,17 17 IF(IHLF-10)11,18,18 C C NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE. 18 IHLF=11 X=X+H GOTO 4 C C THERE IS SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS. 19 X=X+H CALL FCT(X,Y,DERY) DO 20 I=1,NDIM AUX(3,I)=Y(I) 20 AUX(10,I)=DERY(I) N=3 ISW=4 GOTO 100 C 21 N=1 X=X+H CALL FCT(X,Y,DERY) X=PRMT(1) DO 22 I=1,NDIM AUX(11,I)=DERY(I) 220Y(I)=AUX(1,I)+H*(.375D0*AUX(8,I)+.7916666666666667D0*AUX(9,I) 1-.20833333333333333D0*AUX(10,I)+.041666666666666667D0*DERY(I)) 23 X=X+H N=N+1 CALL FCT(X,Y,DERY) CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))6,24,6 24 IF(N-4)25,200,200 25 DO 26 I=1,NDIM AUX(N,I)=Y(I) 26 AUX(N+7,I)=DERY(I) IF(N-3)27,29,200 C 27 DO 28 I=1,NDIM DELT=AUX(9,I)+AUX(9,I) DELT=DELT+DELT 28 Y(I)=AUX(1,I)+.33333333333333333D0*H*(AUX(8,I)+DELT+AUX(10,I)) GOTO 23 C 29 DO 30 I=1,NDIM DELT=AUX(9,I)+AUX(10,I) DELT=DELT+DELT+DELT 30 Y(I)=AUX(1,I)+.375D0*H*(AUX(8,I)+DELT+AUX(11,I)) GOTO 23 C C THE FOLLOWING PART OF SUBROUTINE DHPCG COMPUTES BY MEANS OF C RUNGE-KUTTA METHOD STARTING VALUES FOR THE NOT SELF-STARTING C PREDICTOR-CORRECTOR METHOD. 100 DO 101 I=1,NDIM Z=H*AUX(N+7,I) AUX(5,I)=Z 101 Y(I)=AUX(N,I)+.4D0*Z C Z IS AN AUXILIARY STORAGE LOCATION C Z=X+.4D0*H CALL FCT(Z,Y,DERY) DO 102 I=1,NDIM Z=H*DERY(I) AUX(6,I)=Z 102 Y(I)=AUX(N,I)+.29697760924775360D0*AUX(5,I)+.15875964497103583D0*Z C Z=X+.45573725421878943D0*H CALL FCT(Z,Y,DERY) DO 103 I=1,NDIM Z=H*DERY(I) AUX(7,I)=Z 103 Y(I)=AUX(N,I)+.21810038822592047D0*AUX(5,I)-3.0509651486929308D0* 1AUX(6,I)+3.8328647604670103D0*Z C Z=X+H CALL FCT(Z,Y,DERY) DO 104 I=1,NDIM 1040Y(I)=AUX(N,I)+.17476028226269037D0*AUX(5,I)-.55148066287873294D0* 1AUX(6,I)+1.2055355993965235D0*AUX(7,I)+.17118478121951903D0* 2H*DERY(I) GOTO(9,13,15,21),ISW C C POSSIBLE BREAK-POINT FOR LINKAGE C C STARTING VALUES ARE COMPUTED. C NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD. 200 ISTEP=3 IPAR=ISTEP 201 IF(N-8)204,202,204 C C N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS 202 DO 203 N=2,7 DO 203 I=1,NDIM AUX(N-1,I)=AUX(N,I) 203 AUX(N+6,I)=AUX(N+7,I) N=7 C C N LESS THAN 8 CAUSES N+1 TO GET N 204 N=N+1 C C COMPUTATION OF NEXT VECTOR Y DO 205 I=1,NDIM AUX(N-1,I)=Y(I) 205 AUX(N+6,I)=DERY(I) X=X+H 206 ISTEP=ISTEP+1 IPAR=IPAR+1 DO 207 I=1,NDIM 0DELT=AUX(N-4,I)+1.3333333333333333D0*H*(AUX(N+6,I)+AUX(N+6,I)- 1AUX(N+5,I)+AUX(N+4,I)+AUX(N+4,I)) Y(I)=DELT-.9256198347107438D0*AUX(16,I) 207 AUX(16,I)=DELT C PREDICTOR IS NOW GENERATED IN ROW 16 OF AUX, MODIFIED PREDICTOR C IS GENERATED IN Y. DELT MEANS AN AUXILIARY STORAGE. C CALL FCT(X,Y,DERY) C DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY C DO 208 I=1,NDIM 0DELT=.125D0*(9.D0*AUX(N-1,I)-AUX(N-3,I)+3.D0*H*(DERY(I)+AUX(N+6,I) 1+AUX(N+6,I)-AUX(N+5,I))) AUX(16,I)=AUX(16,I)-DELT 208 Y(I)=DELT+.07438016528925620D0*AUX(16,I) C C TEST WHETHER H MUST BE HALVED OR DOUBLED IF ((X-PRMT(2))/H.GT.1.D-6) GOTO 222 DELT=0.D0 DO 209 I=1,NDIM 209 DELT=DELT+AUX(15,I)*DABS(AUX(16,I)) IF(DELT-PRMT(4))210,222,222 C C H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD. 210 CALL FCT(X,Y,DERY) CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))212,211,212 211 IF(IHLF-11)213,212,212 212 RETURN 213 IF(H*(X-PRMT(2)))214,212,212 214 IF(DABS(X-PRMT(2))-.1D0*DABS(H))212,215,215 215 IF(DELT-.02D0*PRMT(4))217,217,201 C C C H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE C AVAILABLE 217 IF(N-7)201,218,218 218 IF(ISTEP-4)201,219,219 219 IMOD=IPAR/2 IF(IPAR-IMOD-IMOD)201,220,201 220 H=H+H IHLF=IHLF-1 ISTEP=0 IPAR=IMOD DO 221 I=1,NDIM AUX(N-1,I)=AUX(N-2,I) AUX(N-2,I)=AUX(N-4,I) AUX(N-3,I)=AUX(N-6,I) AUX(N+6,I)=AUX(N+5,I) AUX(N+5,I)=AUX(N+3,I) AUX(N+4,I)=AUX(N+1,I) DELT=AUX(N+6,I)+AUX(N+5,I) DELT=DELT+DELT+DELT 2210AUX(16,I)=8.962962962962963D0*(Y(I)-AUX(N-3,I)) 1-3.3611111111111111D0*H*(DERY(I)+DELT+AUX(N+4,I)) GOTO 201 C C C H MUST BE HALVED 222 IHLF=IHLF+1 IF(IHLF-10)223,223,210 223 H=.5D0*H ISTEP=0 IPAR=IPAR+IPAR-2 DO 224 I=1,NDIM 0Y(I)=.390625D-2*(8.D1*AUX(N-1,I)+135.D0*AUX(N-2,I)+4.D1*AUX(N-3,I) 1+AUX(N-4,I))-.1171875D0*(AUX(N+6,I)-6.D0*AUX(N+5,I)-AUX(N+4,I))*H 0AUX(N-4,I)=.390625D-2*(12.D0*AUX(N-1,I)+135.D0*AUX(N-2,I)+ 1108.D0*AUX(N-3,I)+AUX(N-4,I))-.0234375D0*(AUX(N+6,I)+ 218.D0*AUX(N+5,I)-9.D0*AUX(N+4,I))*H AUX(N-3,I)=AUX(N-2,I) 224 AUX(N+4,I)=AUX(N+5,I) X=X-H DELT=X-(H+H) CALL FCT(DELT,Y,DERY) DO 225 I=1,NDIM AUX(N-2,I)=Y(I) AUX(N+5,I)=DERY(I) 225 Y(I)=AUX(N-4,I) DELT=DELT-(H+H) CALL FCT(DELT,Y,DERY) DO 226 I=1,NDIM DELT=AUX(N+5,I)+AUX(N+4,I) DELT=DELT+DELT+DELT 0AUX(16,I)=8.962962962962963D0*(AUX(N-1,I)-Y(I)) 1-3.3611111111111111D0*H*(AUX(N+6,I)+DELT+DERY(I)) 226 AUX(N+3,I)=DERY(I) GOTO 206 END