Trailing-Edge
-
PDP-10 Archives
-
decus_20tap5_198111
-
decus/20-0137/difeq/difeq2.for
There are 2 other files named difeq2.for in the archive. Click here to see a list.
C WESTERN MICHIGAN UNIVERSITY
C DIFEQ2.FOR (FILE NAME ON LIBRARY DECTAPE)
C DIFEQ2.FOR IS CALLED (BY RUNUUO) FROM DIFEQ1.FOR
C FORWMU PROGS. USED: DELETE, RUNUUO
C INTERNAL SUBR. USED: OUTP, HPCG, RKGS
C EXTERNAL SUBR. USED: (GENERATED BY DIFEQ1) FCT.F4, DIV.DAT
C FCT CONTAINS USER SUPPLIED EQS. IN FORTRAN FORM.
C DIV.DAT CONTAINS USER RESPONSES TO DIALOGUE, I.E. NO. OF EQS.
C INTERVAL OF SOL., INITIAL INCREMENT OF X, UPPER ERROR BOUND,
C INITIAL VALUES.
C ABOVE COMMENTS AND RIGHT ADJUSTED COMMENTS PUT IN BY WG
C
DIMENSION KL(8),O(302),N1(10),N2(10)
DIMENSION PRMT(5), Y(20), DERY(20), AUX(26,20)
EXTERNAL FCT
EXTERNAL OUTP
INTEGER O
COMMON /OUTNUM/JDOUT
C ********THIS DATA STATEMENT MUST BE CHANGED IF PROGRAM
C IS NOT EXECUTED FROM 1,4
C FOR SYSTEM OPERATION IT READS: RUN DIFEQ[1,4]
DATA KL/'R ','DIFEQ',0,0,0,0,0,0/
DATA NPR/5H' /
DATA N1/' ',' ',' ',' ',' ',' ',' ',' ',' ','1'/
DATA N2/'1','2','3','4','5','6','7','8','9','0'/
GO TO 301
108 FORMAT(2I)
100 TYPE 102
102 FORMAT(' ENTER INTERVAL OF SOLUTION'/)
ACCEPT 103,PRMT(1),PRMT(2)
103 FORMAT(2F)
TYPE 104
104 FORMAT(' ENTER INITIAL INCREMENT OF X'/)
ACCEPT 103,PRMT(3)
TYPE 105
105 FORMAT(' ENTER UPPER ERROR BOUND'/)
ACCEPT 103,PRMT(4)
TYPE 106
106 FORMAT(' ENTER INITIAL VALUES'/)
ACCEPT 109,(Y(I),I=1,NDIM)
109 FORMAT(10F)
CALL RELEAS(20)
CALL OFILE(20,'DIF')
WRITE(20)PRMT,NDIM,Y
GO TO 302
CALL RELEAS(20)
301 CALL IFILE(20,'DIF')
READ(20)PRMT,NDIM,Y
302 DO 110 I=1,NDIM
110 DERY(I)=1.0/FLOAT(NDIM)
107 TYPE 101
101 FORMAT(' ENTER METHOD OF SOLUTION:',/,' 1 FOR HAMMINGS PREDICT',
1'OR-CORRECTOR METHOD,',/,' 2 FOR 4TH ORDER RUNGE-KUTTA METHOD',/)
ACCEPT 108,MTH,IOUT
IF(MTH.LT.1.OR.MTH.GT.2)GO TO 107
C OUTPUT TO DISK IF IOUT NOT EQUAL ZERO.
JDOUT=21
IF(IOUT.NE.0)GO TO 400
JDOUT=30
WRITE(30,200)
200 FORMAT(///' SOLUTION:'//)
NN=28*NDIM+14
DO 201 I=1,NN
201 O(I)=' '
O(8)='X'
DO 202 I=1,NDIM
K=14*I
O(K+7)='Y'
O(K+8)='('
O(K+9)=N1(I)
O(K+10)=N2(I)
O(K+11)=')'
K=NDIM*14+14*I
O(K+6)='Y'
O(K+7)=NPR
O(K+8)='('
O(K+9)=N1(I)
O(K+10)=N2(I)
202 O(K+11)=')'
WRITE(30,203)(O(J),J=1,NN)
203 FORMAT(70A1)
400 IF(MTH.EQ.1)CALL HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
IF(MTH.EQ.2)CALL RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
IF(IOUT.NE.0)WRITE(JDOUT,402)
402 FORMAT(' -9999999.')
IF(IOUT.NE.0)TYPE 404
404 FORMAT(' OUTPUT FILE IS "FOR21.',
1 'DAT". FORMAT IS "(2X,5G14.7)".',/)
IF(IHLF.EQ.11)TYPE 120
IF(IHLF.GT.11)TYPE 121
TYPE 122
122 FORMAT(//,' TYPE 1 TO ENTER NEW EQUATIONS, 2 TO RECOMPUTE, ',/,
1' 3 TO CHANGE METHOD OF SOLUTION, 4 TO EXIT',/)
ACCEPT 108,KTO
GO TO (300,100,301,126),KTO
126 CALL DELETE('FCT.F4 ')
CALL DELETE ('DIF.DAT ')
CALL EXIT
120 FORMAT(' MORE THAN 10 BISECTIONS OF INCREMENT'/)
121 FORMAT(' INITIAL INCREMENT = 0 OR HAS THE WRONG SIGN'/)
C---------------SEE DATA KL ABOVE. DIFEQ1.EXE IS IN SYS:
300 CALL RUNUUO(KL)
END
C
C SUBROUTINE OUTP
C
C USER SUPPLIED SUBROUTINE TO PRINT SOLUTIN TABLE
C
C---------------Y, DERY, NDIM, X ARE INPUT. JDOUT IS INPUT THRU
C--------------- COMMON /OUTNUM/. IHLF, PRMT APPARENTLY ARE NOT USED.
SUBROUTINE OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
COMMON/OUTNUM/JDOUT
DIMENSION Y(1),DERY(1),PRMT(1)
WRITE(JDOUT,100)X,(Y(I),I=1,NDIM),(DERY(I),I=1,NDIM)
100 FORMAT(2X,5G14.7)
RETURN
END
C
C
C ..................................................................
C SUBROUTINE HPCG
C
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
C---------------PRMT, Y, DERY, NDIM, FCT, OUTP ARE INPUT.
C--------------- IHLF, AUX, RETURNED. PRMT IS MODIFIED.
SUBROUTINE HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
C
C
DIMENSION PRMT(1),Y(1),DERY(1),AUX(16,1)
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
C SUBROUTINE RKGS
C
C PURPOSE
C TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL
C EQUATIONS WITH GIVEN INITIAL VALUES.
C
C USAGE
C CALL RKGS (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 RKGS. EXCEPT PRMT(5)
C THE COMPONENTS ARE NOT DESTROYED BY SUBROUTINE
C RKGS 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 RKGS INITIALIZES
C PRMT(5)=0. IF THE USER WANTS TO TERMINATE
C SUBROUTINE RKGS 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 RKGS DOES NOT REQUIRE
C AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL
C FOR HANDING RESULT VALUES TO THE MAIN PROGRAM
C (CALLING RKGS) 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 RKGS RETURNS WITH
C ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM. ERROR
C 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. THIS
C SUBROUTINE COMPUTES THE RIGHT HAND SIDES DERY OF
C THE SYSTEM TO GIVEN VALUES X AND Y. ITS PARAMETER
C LIST MUST BE X,Y,DERY. SUBROUTINE FCT SHOULD
C NOT 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 RKGS IS TERMINATED.
C AUX - AN AUXILIARY STORAGE ARRAY WITH 8 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 FOURTH ORDER RUNGE-KUTTA
C FORMULAE IN THE MODIFICATION DUE TO GILL. ACCURACY IS
C TESTED COMPARING THE RESULTS OF THE PROCEDURE WITH SINGLE
C AND DOUBLE INCREMENT.
C SUBROUTINE RKGS AUTOMATICALLY ADJUSTS THE INCREMENT DURING
C THE WHOLE COMPUTATION BY HALVING OR DOUBLING. IF MORE THAN
C 10 BISECTIONS OF THE INCREMENT ARE NECESSARY TO GET
C SATISFACTORY ACCURACY, THE SUBROUTINE RETURNS WITH
C ERROR MESSAGE IHLF=11 INTO MAIN PROGRAM.
C TO GET FULL FLEXIBILITY IN OUTPUT, AN OUTPUT SUBROUTINE
C MUST BE FURNISHED BY THE USER.
C FOR REFERENCE, SEE
C RALSTON/WILF, MATHEMATICAL METHODS FOR DIGITAL COMPUTERS,
C WILEY, NEW YORK/LONDON, 1960, PP.110-120.
C
C ..................................................................
C
C---------------PRMT, Y, DERY, NDIM, FCT, OUPT ARE INPUT. IHLF,
C--------------- AUX ARE OUTPUT. PRMT IS MODIFIED.
SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
C
C
DIMENSION Y(1),DERY(1),AUX(8,1),A(4),B(4),C(4),PRMT(1)
DO 1 I=1,NDIM
1 AUX(8,I)=.06666667*DERY(I)
X=PRMT(1)
XEND=PRMT(2)
H=PRMT(3)
PRMT(5)=0.
CALL FCT(X,Y,DERY)
C
C ERROR TEST
IF(H*(XEND-X))38,37,2
C
C PREPARATIONS FOR RUNGE-KUTTA METHOD
2 A(1)=.5
A(2)=.2928932
A(3)=1.707107
A(4)=.1666667
B(1)=2.
B(2)=1.
B(3)=1.
B(4)=2.
C(1)=.5
C(2)=.2928932
C(3)=1.707107
C(4)=.5
C
C PREPARATIONS OF FIRST RUNGE-KUTTA STEP
DO 3 I=1,NDIM
AUX(1,I)=Y(I)
AUX(2,I)=DERY(I)
AUX(3,I)=0.
3 AUX(6,I)=0.
IREC=0
H=H+H
IHLF=-1
ISTEP=0
IEND=0
C
C
C START OF A RUNGE-KUTTA STEP
4 IF((X+H-XEND)*H)7,6,5
5 H=XEND-X
6 IEND=1
C
C RECORDING OF INITIAL VALUES OF THIS STEP
7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT)
IF(PRMT(5))40,8,40
8 ITEST=0
9 ISTEP=ISTEP+1
C
C
C START OF INNERMOST RUNGE-KUTTA LOOP
J=1
10 AJ=A(J)
BJ=B(J)
CJ=C(J)
DO 11 I=1,NDIM
R1=H*DERY(I)
R2=AJ*(R1-BJ*AUX(6,I))
Y(I)=Y(I)+R2
R2=R2+R2+R2
11 AUX(6,I)=AUX(6,I)+R2-CJ*R1
IF(J-4)12,15,15
12 J=J+1
IF(J-3)13,14,13
13 X=X+.5*H
14 CALL FCT(X,Y,DERY)
GOTO 10
C END OF INNERMOST RUNGE-KUTTA LOOP
C
C
C TEST OF ACCURACY
15 IF(ITEST)16,16,20
C
C IN CASE ITEST=0 THERE IS NO POSSIBILITY FOR TESTING OF ACCURACY
16 DO 17 I=1,NDIM
17 AUX(4,I)=Y(I)
ITEST=1
ISTEP=ISTEP+ISTEP-2
18 IHLF=IHLF+1
X=X-H
H=.5*H
DO 19 I=1,NDIM
Y(I)=AUX(1,I)
DERY(I)=AUX(2,I)
19 AUX(6,I)=AUX(3,I)
GOTO 9
C
C IN CASE ITEST=1 TESTING OF ACCURACY IS POSSIBLE
20 IMOD=ISTEP/2
IF(ISTEP-IMOD-IMOD)21,23,21
21 CALL FCT(X,Y,DERY)
DO 22 I=1,NDIM
AUX(5,I)=Y(I)
22 AUX(7,I)=DERY(I)
GOTO 9
C
C COMPUTATION OF TEST VALUE DELT
23 DELT=0.
DO 24 I=1,NDIM
24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I))
IF(DELT-PRMT(4))28,28,25
C
C ERROR IS TOO GREAT
25 IF(IHLF-10)26,36,36
26 DO 27 I=1,NDIM
27 AUX(4,I)=AUX(5,I)
ISTEP=ISTEP+ISTEP-4
X=X-H
IEND=0
GOTO 18
C
C RESULT VALUES ARE GOOD
28 CALL FCT(X,Y,DERY)
DO 29 I=1,NDIM
AUX(1,I)=Y(I)
AUX(2,I)=DERY(I)
AUX(3,I)=AUX(6,I)
Y(I)=AUX(5,I)
29 DERY(I)=AUX(7,I)
CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT)
IF(PRMT(5))40,30,40
30 DO 31 I=1,NDIM
Y(I)=AUX(1,I)
31 DERY(I)=AUX(2,I)
IREC=IHLF
IF(IEND)32,32,39
C
C INCREMENT GETS DOUBLED
32 IHLF=IHLF-1
ISTEP=ISTEP/2
H=H+H
IF(IHLF)4,33,33
33 IMOD=ISTEP/2
IF(ISTEP-IMOD-IMOD)4,34,4
34 IF(DELT-.02*PRMT(4))35,35,4
35 IHLF=IHLF-1
ISTEP=ISTEP/2
H=H+H
GOTO 4
C
C
C RETURNS TO CALLING PROGRAM
36 IHLF=11
CALL FCT(X,Y,DERY)
GOTO 39
37 IHLF=12
GOTO 39
38 IHLF=13
39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT)
40 RETURN
END