Directory of image this file is from
This file as a plain text file
C
C ..................................................................
C
C SUBROUTINE DRKGS
C
C PURPOSE
C TO SOLVE A SYSTEM OF FIRST ORDER ORDINARY DIFFERENTIAL
C EQUATIONS WITH GIVEN INITIAL VALUES.
C
C USAGE
C CALL DRKGS (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 DRKGS. EXCEPT PRMT(5) THE COMPONENTS
C ARE NOT DESTROYED BY SUBROUTINE DRKGS 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 DRKGS INITIALIZES
C PRMT(5)=0. IF THE USER WANTS TO TERMINATE
C SUBROUTINE DRKGS 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 DRKGS DOES NOT REQUIRE
C AND CHANGE THEM. NEVERTHELESS THEY MAY BE USEFUL
C FOR HANDING RESULT VALUES TO THE MAIN PROGRAM
C (CALLING DRKGS) 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 DRKGS 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 DRKGS IS TERMINATED.
C AUX - DOUBLE PRECISION AUXILIARY STORAGE ARRAY WITH 8
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 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 DRKGS 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
SUBROUTINE DRKGS(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)
DOUBLE PRECISION PRMT,Y,DERY,AUX,A,B,C,X,XEND,H,AJ,BJ,CJ,R1,R2,
1DELT
DO 1 I=1,NDIM
1 AUX(8,I)=.066666666666666667D0*DERY(I)
X=PRMT(1)
XEND=PRMT(2)
H=PRMT(3)
PRMT(5)=0.D0
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)=.5D0
A(2)=.29289321881345248D0
A(3)=1.7071067811865475D0
A(4)=.16666666666666667D0
B(1)=2.D0
B(2)=1.D0
B(3)=1.D0
B(4)=2.D0
C(1)=.5D0
C(2)=.29289321881345248D0
C(3)=1.7071067811865475D0
C(4)=.5D0
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.D0
3 AUX(6,I)=0.D0
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+.5D0*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=.5D0*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.D0
DO 24 I=1,NDIM
24 DELT=DELT+AUX(8,I)*DABS(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-.02D0*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