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