C #RKSTB V1A 19-JAN-73 C LAST UPDATE: C SUBROUTINE RKSTB(DERIV,Y0,YP0,Y1,N,X,H,IX,WK) DIMENSION Y0(N), YP0(N), Y1(N), WK(N,3) C C *PURPOSE. C TO MAKE A SINGLE STEP IN INTEGRATING A SET OF ORDINARY DIFFERENTIAL C EQUATIONS USING THE FOURTH ORDER STANDARD RUNGE-KUTTA METHOD. C C *PARAMETERS: C C DERIV- IS THE EXTERNAL NAME OF THE USER'S ROUTINE TO EVALUATE THE C RIGHT HAND SIDES OF THE DIFFERENTIAL EQUATIONS. CALLED AS C CALL DERIV(Y,YP,N,X) C WHERE: C Y - IS A REAL ARRAY OF LENGTH "N" CONTAINING THE DEPENDENT C VARIABLES Y(I). C YP - IS A REAL ARRAY IN WHICH THE DERIVATIVES OF THE Y(I) C ARE RETURNED BY "DERIV" C N - IS THE NUMBER OF EQUATIONS BEING INTEGRATED. C X - IS THE INDEPENDENT VARIABLE. C Y0 - IS THE INITIAL VALUE FOR Y(I). A REAL ARRAY OF LENGTH N C YP0 - IS THE INITIAL VALUE OF THE DERIVATIVES OF Y(I) AT X, C INPUT IF IX.NE.0 AND OUTPUT OTHERWISE. C Y1 - IS THE ARRAY OF FINAL Y(I) VALUES AT X+H C N - IS THE NUMBER OF EQUATIONS BEING INTEGRATED. C X - IS THE INDEPENDENT VARIABLE RETURNED AS X+H C H - IS THE FIXED SINGLE STEP TAKEN IN X. C IX - IS ZERO IF THE DERIVATIVES AT THE INPUT VALUE X HAVE TO C BE CALCULATED, ELSE IS NON-ZERO AND THE DERIVATIVES ARE C INPUT IN YP0. C WK - IS A WORKING ARRAY OF AT LEAST 3*N REAL ELEMENTS. C C INPUT VARIABLES: DERIV, Y0, N, X, H, IX AND POSSIBLY YP0 C OUTPUT VARIABLES: YP0, Y1, X C C NOTE. THE VECTORS Y, YP0 AND Y1 MUST ALL BE DISTINCT - AND DISTINCT C FROM THE COLUMNS OF WK!! C C *METHOD. C THE FOURTH ORDER RUNGE-KUTTA ALGORITHM IS CARRIED IN THE C STRAIGHTFORWARD WAY. THE CODING IS DESIGNED TO REDUCE TO A MINIMUM C THE NUMBER OF FLOATING POINT MULTIPLICATIONS. C IF IX IS ZERO THE DERIVATIVES AT THE INITIAL X VALUE ARE CALCULATED C IN YP0. OTHERWISE THEY ARE ASSUMED TO BE THERE ALREADY. C C *ACCURACY. C ROUNDING ERROR WILL BE OF THE ORDER OF THE MACHINE SIGNIFICANCE C TRUNCATION ERRORS DEPEND ON THE STEP SIZE AND THE PROBLEM. C C *RESTRICTIONS. C C *ERROR CONDITIONS. C C *NON STANDARD ROUTINES CALLED. C DERIV(USER'S) C C *TYPICAL TIMES. C C *ORIGIN. M.R.MANNING. C C *COMMENTS. C C #END. C C C #START. HBY2=H/2. HBY6=H/6. C IF(IX) 10,5,10 5 CALL DERIV(Y0,YP0,N,X) C 10 DO 12 I=1,N 12 Y1(I)= Y0(I) + HBY2*YP0(I) XA=X+HBY2 CALL DERIV(Y1,WK,N,XA) DO 14 I=1,N 14 Y1(I)= Y0(I) + HBY2*WK(I,1) CALL DERIV(Y1,WK(1,2), N,XA) DO 16 I=1,N 16 Y1(I)= Y0(I) + H*WK(I,2) X=X+H CALL DERIV(Y1,WK(1,3), N,X) DO 20 I=1,N T=WK(I,1) + WK(I,2) T=HBY6*(YP0(I) + T + T +WK(I,3) ) WK(I,1)=T 20 Y1(I)=Y0(I) + T C RETURN END