File: RK2.FT of Tape: Various/ETH/eth11-1

(Source file text) 

---

C     ..................................................................
C
C        SUBROUTINE RK2
C
C        PURPOSE
C           INTEGRATES A FIRST ORDER DIFFERENTIAL EQUATION
C           DY/DX=FUN(X,Y) AND PRODUCES A TABLE OF INTEGRATED VALUES
C
C        USAGE
C           CALL RK2(FUN,H,XI,YI,K,N,VEC)
C
C        DESCRIPTION OF PARAMETERS
C           FUN-USER-SUPPLIED FUNCTION SUBPROGRAM WITH ARGUMENTS X,Y
C               WHICH GIVES DY/DX
C           H  -STEP SIZE
C           XI -INITIAL VALUE OF X
C           YI -INITIAL VALUE OF Y WHERE YI=Y(XI)
C           K  -THE INTERVAL AT WHICH COMPUTED VALUES ARE TO BE STORED
C           N  -THE NUMBER OF VALUES TO BE STORED
C           VEC-THE RESULTANT VECTOR OF LENGTH N IN WHICH COMPUTED
C               VALUES OF Y ARE TO BE STORED
C
C        REMARKS
C           NONE
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           FUN - USER-SUPPLIED FUNCTION SUBPROGRAM FOR DY/DX
C           CALLING PROGRAM MUST HAVE FORTRAN EXTERNAL STATEMENT
C           CONTAINING NAMES OF FUNCTION SUBPROGRAMS LISTED IN CALL TO
C           RK2
C
C        METHOD
C           FOURTH ORDER RUNGE-KUTTA INTEGRATION ON A RECURSIVE BASIS AS
C           SHOWN IN F.B. HILDEBRAND, 'INTRODUCTION TO NUMERICAL
C           ANALYSIS', MCGRAW-HILL, NEW YORK, 1956
C
C     ..................................................................
C
      SUBROUTINE RK2(FUN,H,XI,YI,K,N,VEC)
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION H,XI,YI,VEC,H2,Y,X,T1,T2,T3,T4,FUN
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        USER FUNCTION SUBPROGRAM, FUN, MUST BE IN DOUBLE PRECISION.
C
C        ...............................................................
C
      DIMENSION VEC(1)
      H2=H/2.
      Y=YI
      X=XI
      DO 2 I=1,N
      DO 1 J=1,K
      T1=H*FUN(X,Y)
      T2=H*FUN(X+H2,Y+T1/2.)
      T3=H*FUN(X+H2,Y+T2/2.)
      T4=H*FUN(X+H,Y+T3)
      Y= Y+(T1+2.*T2+2.*T3+T4)/6.
    1 X=X+H
    2 VEC(I)=Y
      RETURN
      END
C

---