C SUBROUTINE 'RKGS' 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 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 SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) C C DIMENSION Y(*),DERY(*),AUX(8,*),A(4),B(4),C(4),PRMT(*) 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 C C======================================================================= C