Directory of image this file is from
This file as a plain text file
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