File: RK1.FT of Tape: Various/ETH/eth11-1
(Source file text)
C ..................................................................
C
C SUBROUTINE RK1
C
C PURPOSE
C INTEGRATES A FIRST ORDER DIFFERENTIAL EQUATION
C DY/DX=FUN(X,Y) UP TO A SPECIFIED FINAL VALUE
C
C USAGE
C CALL RK1(FUN,HI,XI,YI,XF,YF,ANSX,ANSY,IER)
C
C DESCRIPTION OF PARAMETERS
C FUN -USER-SUPPLIED FUNCTION SUBPROGRAM WITH ARGUMENTS X,Y
C WHICH GIVES DY/DX
C HI -THE STEP SIZE
C XI -INITIAL VALUE OF X
C YI -INITIAL VALUE OF Y WHERE YI=Y(XI)
C XF -FINAL VALUE OF X
C YF -FINAL VALUE OF Y
C ANSX-RESULTANT FINAL VALUE OF X
C ANSY-RESULTANT FINAL VALUE OF Y
C EITHER ANSX WILL EQUAL XF OR ANSY WILL EQUAL YF
C DEPENDING ON WHICH IS REACHED FIRST
C IER -ERROR CODE
C IER=0 NO ERROR
C IER=1 STEP SIZE IS ZERO
C
C REMARKS
C IF XI IS GREATER THAN XF, ANSX=XI AND ANSY=YI
C IF H IS ZERO, IER IS SET TO ONE, ANSX IS SET TO XI, AND
C ANSY IS SET TO ZERO
C
C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C FUN IS A TWO ARGUMENT FUNCTION SUBPROGRAM FURNISHED BY THE
C USER. DY/DX=FUN (X,Y)
C CALLING PROGRAM MUST HAVE FORTRAN EXTERNAL STATEMENT
C CONTAINING NAMES OF FUNCTION SUBPROGRAMS LISTED IN CALL TO
C RK1
C
C METHOD
C USES FOURTH ORDER RUNGE-KUTTA INTEGRATION PROCESS ON A
C RECURSIVE BASIS AS SHOWN IN F.B. HILDEBRAND, 'INTRODUCTION
C TO NUMERICAL ANALYSIS',MCGRAW-HILL,1956. PROCESS IS
C TERMINATED AND FINAL VALUE ADJUSTED WHEN EITHER XF OR YF
C IS REACHED.
C
C ..................................................................
C
SUBROUTINE RK1(FUN,HI,XI,YI,XF,YF,ANSX,ANSY,IER)
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 HI,XI,YI,XF,YF,ANSX,ANSY,H,XN,YN,HNEW,XN1,YN1,
C 1 XX,YY,XNEW,YNEW,H2,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
C IF XF IS LESS THAN OR EQUAL TO XI, RETURN XI,YI AS ANSWER
C
IER=0
IF(XF-XI) 11,11,12
11 ANSX=XI
ANSY=YI
RETURN
C
C TEST INTERVAL VALUE
C
12 H=HI
IF(HI) 16,14,20
14 IER=1
ANSX=XI
ANSY=0.0
RETURN
16 H=-HI
C
C SET XN=INITIAL X,YN=INITIAL Y
C
20 XN=XI
YN=YI
C
C INTEGRATE ONE TIME STEP
C
HNEW=H
JUMP=1
GO TO 170
25 XN1=XX
YN1=YY
C
C COMPARE XN1 (=X(N+1)) TO X FINAL AND BRANCH ACCORDINGLY
C
IF(XN1-XF)50,30,40
C
C XN1=XF, RETURN (XF,YN1) AS ANSWER
C
30 ANSX=XF
ANSY=YN1
GO TO 160
C
C XN1 GREATER THAN XF, SET NEW STEP SIZE AND INTEGRATE ONE STEP
C RETURN RESULTS OF INTEGRATION AS ANSWER
C
40 HNEW=XF-XN
JUMP=2
GO TO 170
45 ANSX=XX
ANSY=YY
GO TO 160
C
C XN1 LESS THAN X FINAL, CHECK IF (YN,YN1) SPAN Y FINAL
C
C
50 IF((YN1-YF)*(YF-YN))60,70,110
C
C YN1 AND YN DO NOT SPAN YF. SET (XN,YN) AS (XN1,YN1) AND REPEAT
C
60 YN=YN1
XN=XN1
GO TO 170
C
C EITHER YN OR YN1 =YF. CHECK WHICH AND SET PROPER (X,Y) AS ANSWER
C
70 IF(YN1-YF)80,100,80
80 ANSY=YN
ANSX=XN
GO TO 160
100 ANSY=YN1
ANSX=XN1
GO TO 160
C
C YN AND YN1 SPAN YF. TRY TO FIND X VALUE ASSOCIATED WITH YF
C
110 DO 140 I=1,10
C
C INTERPOLATE TO FIND NEW TIME STEP AND INTEGRATE ONE STEP
C TRY TEN INTERPOLATIONS AT MOST
C
HNEW=((YF-YN )/(YN1-YN))*(XN1-XN)
JUMP=3
GO TO 170
115 XNEW=XX
YNEW=YY
C
C COMPARE COMPUTED Y VALUE WITH YF AND BRANCH
C
IF(YNEW-YF)120,150,130
C
C ADVANCE, YF IS BETWEEN YNEW AND YN1
C
120 YN=YNEW
XN=XNEW
GO TO 140
C
C ADVANCE, YF IS BETWEEN YN AND YNEW
C
130 YN1=YNEW
XN1=XNEW
140 CONTINUE
C
C RETURN (XNEW,YF) AS ANSWER
C
150 ANSX=XNEW
ANSY=YF
160 RETURN
C
170 H2=HNEW/2.0
T1=HNEW*FUN(XN,YN)
T2=HNEW*FUN(XN+H2,YN+T1/2.0)
T3=HNEW*FUN(XN+H2,YN+T2/2.0)
T4=HNEW*FUN(XN+HNEW,YN+T3)
YY=YN+(T1+2.0*T2+2.0*T3+T4)/6.0
XX=XN+HNEW
GO TO (25,45,115), JUMP
C
END
C