C SUBROUTINE 'HPCG' FROM THE IBM SCIENTIFIC SUBROUTINE PACKAGE. C C NOTE: TO CONFORM WITH THE FORTRAN77 STANDARD, DUMMY ARRAY DIMENSIONS C (1) HAVE BEEN CHANGED TO (*). C C .................................................................. C C SUBROUTINE HPCG C C PURPOSE C TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY GENERAL C DIFFERENTIAL EQUATIONS WITH GIVEN INITIAL VALUES. C C USAGE C CALL HPCG (PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) C PARAMETERS FCT AND OUTP REQUIRE AN EXTERNAL STATEMENT. C C DESCRIPTION OF PARAMETERS C PRMT - AN INPUT AND OUTPUT VECTOR WITH DIMENSION GREATER C OR EQUAL TO 5, WHICH SPECIFIES THE PARAMETERS OF C THE INTERVAL AND OF ACCURACY AND WHICH SERVES FOR C COMMUNICATION BETWEEN OUTPUT SUBROUTINE (FURNISHED C BY THE USER) AND SUBROUTINE HPCG. EXCEPT PRMT(5) C THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE C HPCG 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 HPCG INITIALIZES C PRMT(5)=0. IF THE USER WANTS TO TERMINATE C SUBROUTINE HPCG 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 HPCG DOES NOT REQUIRE C AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL C FOR HANDING RESULT VALUES TO THE MAIN PROGRAM C (CALLING HPCG) WHICH ARE OBTAINED BY SPECIAL C MANIPULATIONS WITH OUTPUT DATA IN SUBROUTINE OUTP. C Y - INPUT VECTOR OF INITIAL VALUES. (DESTROYED) C LATERON Y IS THE RESULTING VECTOR OF DEPENDENT C VARIABLES COMPUTED AT INTERMEDIATE POINTS X. C DERY - INPUT VECTOR OF ERROR WEIGHTS. (DESTROYED) C THE SUM OF ITS COMPONENTS MUST BE EQUAL TO 1. C LATERON DERY IS THE VECTOR OF DERIVATIVES, WHICH C BELONG TO FUNCTION VALUES Y AT A POINT 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 HPCG 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 HPCG IS TERMINATED. C AUX - AN AUXILIARY STORAGE ARRAY WITH 16 ROWS AND NDIM C 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 HPCG 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 HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) C DIMENSION PRMT(*),Y(*),DERY(*),AUX(16,*) N=1 IHLF=0 X=PRMT(1) H=PRMT(3) PRMT(5)=0. DO 1 I=1,NDIM AUX(16,I)=0. 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=.5*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. DO 16 I=1,NDIM 16 DELT=DELT+AUX(15,I)*ABS(Y(I)-AUX(4,I)) DELT=.06666667*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*(.375*AUX(8,I)+.7916667*AUX(9,I) 1-.2083333*AUX(10,I)+.04166667*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)+.3333333*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)+.375*H*(AUX(8,I)+DELT+AUX(11,I)) GOTO 23 C C THE FOLLOWING PART OF SUBROUTINE HPCG 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)+.4*Z C Z IS AN AUXILIARY STORAGE LOCATION C Z=X+.4*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)+.2969776*AUX(5,I)+.1587596*Z C Z=X+.4557372*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)+.2181004*AUX(5,I)-3.050965*AUX(6,I)+3.832865*Z C Z=X+H CALL FCT(Z,Y,DERY) DO 104 I=1,NDIM 1040Y(I)=AUX(N,I)+.1747603*AUX(5,I)-.5514807*AUX(6,I) 1+1.205536*AUX(7,I)+.1711848*H*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 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 DO 207 I=1,NDIM 0DELT=AUX(N-4,I)+1.333333*H*(AUX(N+6,I)+AUX(N+6,I)-AUX(N+5,I)+ 1AUX(N+4,I)+AUX(N+4,I)) Y(I)=DELT-.9256198*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=.125*(9.*AUX(N-1,I)-AUX(N-3,I)+3.*H*(DERY(I)+AUX(N+6,I)+ 1AUX(N+6,I)-AUX(N+5,I))) AUX(16,I)=AUX(16,I)-DELT 208 Y(I)=DELT+.07438017*AUX(16,I) C C TEST WHETHER H MUST BE HALVED OR DOUBLED DELT=0. DO 209 I=1,NDIM 209 DELT=DELT+AUX(15,I)*ABS(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(ABS(X-PRMT(2))-.1*ABS(H))212,215,215 215 IF(DELT-.02*PRMT(4))216,216,201 C C C H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE C AVAILABLE 216 IF(IHLF)201,201,217 217 IF(N-7)201,218,218 218 IF(ISTEP-4)201,219,219 219 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)201,220,201 220 H=H+H IHLF=IHLF-1 ISTEP=0 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.962963*(Y(I)-AUX(N-3,I))-3.361111*H*(DERY(I)+DELT 1+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=.5*H ISTEP=0 DO 224 I=1,NDIM 0Y(I)=.00390625*(80.*AUX(N-1,I)+135.*AUX(N-2,I)+40.*AUX(N-3,I)+ 1AUX(N-4,I))-.1171875*(AUX(N+6,I)-6.*AUX(N+5,I)-AUX(N+4,I))*H 0AUX(N-4,I)=.00390625*(12.*AUX(N-1,I)+135.*AUX(N-2,I)+ 1108.*AUX(N-3,I)+AUX(N-4,I))-.0234375*(AUX(N+6,I)+18.*AUX(N+5,I)- 29.*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.962963*(AUX(N-1,I)-Y(I))-3.361111*H*(AUX(N+6,I)+DELT 1+DERY(I)) 226 AUX(N+3,I)=DERY(I) GOTO 206 END C C======================================================================= C