Power-Programmierung/CD1.mdf/fortran/library/ssp/orddifeq/drkgs.for

Text File | 1985-11-29 | 9KB | 264 lines

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