C C .................................................................. C C SUBROUTINE HPCL C C PURPOSE C TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY LINEAR C DIFFERENTIAL EQUATIONS WITH GIVEN INITIAL VALUES. C C USAGE C CALL HPCL (PRMT,Y,DERY,NDIM,IHLF,AFCT,FCT,OUTP,AUX,A) C PARAMETERS AFCT,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 HPCL. EXCEPT PRMT(5) C THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE C HPCL 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 HPCL INITIALIZES C PRMT(5)=0. IF THE USER WANTS TO TERMINATE C SUBROUTINE HPCL 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 HPCL DOES NOT REQUIRE C AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL C FOR HANDING RESULT VALUES TO THE MAIN PROGRAM C (CALLING HPCL) 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 HPCL 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 AFCT - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT C COMPUTES MATRIX A (FACTOR OF VECTOR Y ON THE C RIGHT HAND SIDE OF THE SYSTEM) FOR A GIVEN X-VALUE. C ITS PARAMETER LIST MUST BE X,A. THE SUBROUTINE C SHOULD NOT DESTROY X. C FCT - THE NAME OF AN EXTERNAL SUBROUTINE USED. IT C COMPUTES VECTOR F (INHOMOGENEOUS PART OF THE C RIGHT HAND SIDE OF THE SYSTEM) FOR A GIVEN X-VALUE. C ITS PARAMETER LIST MUST BE X,F. THE SUBROUTINE C SHOULD NOT DESTROY X. 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 HPCL IS TERMINATED. C AUX - AN AUXILIARY STORAGE ARRAY WITH 16 ROWS AND NDIM C COLUMNS. C A - AN NDIM BY NDIM MATRIX, WHICH IS USED AS AUXILIARY C STORAGE ARRAY. 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 AFCT(X,A), FCT(X,F) 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 HPCL 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 HPCL(PRMT,Y,DERY,NDIM,IHLF,AFCT,FCT,OUTP,AUX,A) C C C THE FOLLOWING FIRST PART OF SUBROUTINE HPCL (UNTIL FIRST BREAK- C POINT FOR LINKAGE) HAS TO STAY IN CORE DURING THE WHOLE C COMPUTATION C DIMENSION PRMT(1),Y(1),DERY(1),AUX(16,1),A(1) GOTO 100 C C THIS PART OF SUBROUTINE HPCL COMPUTES THE RIGHT HAND SIDE DERY OF C THE GIVEN SYSTEM OF LINEAR DIFFERENTIAL EQUATIONS. 1 CALL AFCT(X,A) CALL FCT(X,DERY) DO 3 M=1,NDIM LL=M-NDIM HS=0. DO 2 L=1,NDIM LL=LL+NDIM 2 HS=HS+A(LL)*Y(L) 3 DERY(M)=HS+DERY(M) GOTO(105,202,204,206,115,122,125,308,312,327,329,128),ISW2 C C POSSIBLE BREAK-POINT FOR LINKAGE C 100 N=1 IHLF=0 X=PRMT(1) H=PRMT(3) PRMT(5)=0. DO 101 I=1,NDIM AUX(16,I)=0. AUX(15,I)=DERY(I) 101 AUX(1,I)=Y(I) IF(H*(PRMT(2)-X))103,102,104 C C ERROR RETURNS 102 IHLF=12 GOTO 104 103 IHLF=13 C C COMPUTATION OF DERY FOR STARTING VALUES 104 ISW2=1 GOTO 1 C C RECORDING OF STARTING VALUES 105 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))107,106,107 106 IF(IHLF)108,108,107 107 RETURN 108 DO 109 I=1,NDIM 109 AUX(8,I)=DERY(I) C C COMPUTATION OF AUX(2,I) ISW1=1 GOTO 200 C 110 X=X+H DO 111 I=1,NDIM 111 AUX(2,I)=Y(I) C C INCREMENT H IS TESTED BY MEANS OF BISECTION 112 IHLF=IHLF+1 X=X-H DO 113 I=1,NDIM 113 AUX(4,I)=AUX(2,I) H=.5*H N=1 ISW1=2 GOTO 200 C 114 X=X+H ISW2=5 GOTO 1 115 N=2 DO 116 I=1,NDIM AUX(2,I)=Y(I) 116 AUX(9,I)=DERY(I) ISW1=3 GOTO 200 C C COMPUTATION OF TEST VALUE DELT 117 DELT=0. DO 118 I=1,NDIM 118 DELT=DELT+AUX(15,I)*ABS(Y(I)-AUX(4,I)) DELT=.06666667*DELT IF(DELT-PRMT(4))121,121,119 119 IF(IHLF-10)112,120,120 C C NO SATISFACTORY ACCURACY AFTER 10 BISECTIONS. ERROR MESSAGE. 120 IHLF=11 X=X+H GOTO 104 C C SATISFACTORY ACCURACY AFTER LESS THAN 11 BISECTIONS 121 X=X+H ISW2=6 GOTO 1 122 DO 123 I=1,NDIM AUX(3,I)=Y(I) 123 AUX(10,I)=DERY(I) N=3 ISW1=4 GOTO 200 C 124 N=1 X=X+H ISW2=7 GOTO 1 125 X=PRMT(1) DO 126 I=1,NDIM AUX(11,I)=DERY(I) 126 Y(I)=AUX(1,I)+H*(.375*AUX(8,I)+.7916667*AUX(9,I) 1-.2083333*AUX(10,I)+.04166667*DERY(I)) 127 X=X+H N=N+1 ISW2=12 GOTO 1 128 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))107,129,107 129 IF(N-4)130,300,300 130 DO 131 I=1,NDIM AUX(N,I)=Y(I) 131 AUX(N+7,I)=DERY(I) IF(N-3)132,134,300 C 132 DO 133 I=1,NDIM DELT=AUX(9,I)+AUX(9,I) DELT=DELT+DELT 133 Y(I)=AUX(1,I)+.3333333*H*(AUX(8,I)+DELT+AUX(10,I)) GOTO 127 C 134 DO 135 I=1,NDIM DELT=AUX(9,I)+AUX(10,I) DELT=DELT+DELT+DELT 135 Y(I)=AUX(1,I)+.375*H*(AUX(8,I)+DELT+AUX(11,I)) GOTO 127 C C THE FOLLOWING PART OF SUBROUTINE HPCL COMPUTES BY MEANS OF C RUNGE-KUTTA METHOD STARTING VALUES FOR THE NOT SELF-STARTING C PREDICTOR-CORRECTOR METHOD. 200 Z=X DO 201 I=1,NDIM X=H*AUX(N+7,I) AUX(5,I)=X 201 Y(I)=AUX(N,I)+.4*X C X IS AN AUXILIARY STORAGE LOCATION C X=Z+.4*H ISW2=2 GOTO 1 202 DO 203 I=1,NDIM X=H*DERY(I) AUX(6,I)=X 203 Y(I)=AUX(N,I)+.2969776*AUX(5,I)+.1587596*X C X=Z+.4557372*H ISW2=3 GOTO 1 204 DO 205 I=1,NDIM X=H*DERY(I) AUX(7,I)=X 205 Y(I)=AUX(N,I)+.2181004*AUX(5,I)-3.050965*AUX(6,I)+3.832865*X C X=Z+H ISW2=4 GOTO 1 206 DO 207 I=1,NDIM 207 Y(I)=AUX(N,I)+.1747603*AUX(5,I)-.5514807*AUX(6,I) 1+1.205536*AUX(7,I)+.1711848*H*DERY(I) X=Z GOTO(110,114,117,124),ISW1 C C POSSIBLE BREAK-POINT FOR LINKAGE C C STARTING VALUES ARE COMPUTED. C NOW START HAMMINGS MODIFIED PREDICTOR-CORRECTOR METHOD. 300 ISTEP=3 301 IF(N-8)304,302,304 C C N=8 CAUSES THE ROWS OF AUX TO CHANGE THEIR STORAGE LOCATIONS 302 DO 303 N=2,7 DO 303 I=1,NDIM AUX(N-1,I)=AUX(N,I) 303 AUX(N+6,I)=AUX(N+7,I) N=7 C C N LESS THAN 8 CAUSES N+1 TO GET N 304 N=N+1 C C COMPUTATION OF NEXT VECTOR Y DO 305 I=1,NDIM AUX(N-1,I)=Y(I) 305 AUX(N+6,I)=DERY(I) X=X+H 306 ISTEP=ISTEP+1 DO 307 I=1,NDIM DELT=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) 307 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. ISW2=8 GOTO 1 C DERIVATIVE OF MODIFIED PREDICTOR IS GENERATED IN DERY C 308 DO 309 I=1,NDIM DELT=.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 309 Y(I)=DELT+.07438017*AUX(16,I) C C TEST WHETHER H MUST BE HALVED OR DOUBLED DELT=0. DO 310 I=1,NDIM 310 DELT=DELT+AUX(15,I)*ABS(AUX(16,I)) IF(DELT-PRMT(4))311,324,324 C C H MUST NOT BE HALVED. THAT MEANS Y(I) ARE GOOD. 311 ISW2=9 GOTO 1 312 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))314,313,314 313 IF(IHLF-11)315,314,314 314 RETURN 315 IF(H*(X-PRMT(2)))316,314,314 316 IF(ABS(X-PRMT(2))-.1*ABS(H))314,317,317 317 IF(DELT-.02*PRMT(4))318,318,301 C C C H COULD BE DOUBLED IF ALL NECESSARY PRECEEDING VALUES ARE C AVAILABLE 318 IF(IHLF)301,301,319 319 IF(N-7)301,320,320 320 IF(ISTEP-4)301,321,321 321 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)301,322,301 322 H=H+H IHLF=IHLF-1 ISTEP=0 DO 323 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 323 AUX(16,I)=8.962963*(Y(I)-AUX(N-3,I))-3.361111*H*(DERY(I)+DELT 1+AUX(N+4,I)) GOTO 301 C C C H MUST BE HALVED 324 IHLF=IHLF+1 IF(IHLF-10)325,325,311 325 H=.5*H ISTEP=0 DO 326 I=1,NDIM Y(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 AUX(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) 326 AUX(N+4,I)=AUX(N+5,I) DELT=X-H X=DELT-(H+H) ISW2=10 GOTO 1 327 DO 328 I=1,NDIM AUX(N-2,I)=Y(I) AUX(N+5,I)=DERY(I) 328 Y(I)=AUX(N-4,I) X=X-(H+H) ISW2=11 GOTO 1 329 X=DELT DO 330 I=1,NDIM DELT=AUX(N+5,I)+AUX(N+4,I) DELT=DELT+DELT+DELT AUX(16,I)=8.962963*(AUX(N-1,I)-Y(I))-3.361111*H*(AUX(N+6,I)+DELT 1+DERY(I)) 330 AUX(N+3,I)=DERY(I) GOTO 306 END