C FIGURE 1.6.1 SUBROUTINE RKF(FSUB,N,T0,TFINAL,HINIT,U,TOL) IMPLICIT DOUBLE PRECISION(A-H,O-Z) C C ARGUMENT DESCRIPTIONS C C FSUB - (INPUT) THE NAME OF A USER-SUPPLIED SUBROUTINE TO EVALUATE C THE RIGHT HAND SIDE(S) OF THE DIFFERENTIAL EQUATIONS C SYSTEM U' = F(T,U). FSUB SHOULD HAVE THE CALLING SEQUENCE C SUBROUTINE FSUB(N,T,U,F) C AND SHOULD EVALUATE F(1)...F(N), GIVEN N,T AND U(1)...U(N). C FSUB MUST BE LISTED IN AN EXTERNAL STATEMENT IN THE C CALLING PROGRAM. C N - (INPUT) THE NUMBER OF DIFFERENTIAL EQUATIONS C T0 - (INPUT) THE INITIAL VALUE OF T. C TFINAL - (INPUT) THE FINAL VALUE OF T. C HINIT - (INPUT) THE INITIAL STEPSIZE TO BE USED. C U - (INPUT) A VECTOR OF LENGTH N, HOLDING THE INITIAL VALUES C OF U AT T=T0. C (OUTPUT) THE SOLUTION VALUES AT T=TFINAL. C THE SOLUTION COMPONENTS ARE PRINTED EACH TIME STEP. C TOL - (INPUT) THE DESIRED BOUND ON THE LOCAL ERROR PER UNIT TIME C OF THE SOLUTION COMPONENTS. C DIMENSION U(N),V1(N),V2(N),V3(N),V4(N),V5(N),V6(N),U4(N),U5(N), & UTEMP(N) TK = T0 H = HINIT 5 TKP1 = MIN(TK+H , TFINAL) H = TKP1-TK C TAKE A STEP WITH STEPSIZE H CALL FSUB(N , TK , U , V1) DO 10 I=1,N 10 UTEMP(I) = U(I) + H*V1(I)/4 CALL FSUB(N , TK+H/4 , UTEMP , V2) DO 15 I=1,N 15 UTEMP(I) = U(I) + H*(3*V1(I)+9*V2(I)) / 32 CALL FSUB(N , TK+3*H/8 , UTEMP , V3) DO 20 I=1,N 20 UTEMP(I) = U(I) + H*(1932*V1(I)-7200*V2(I)+7296*V3(I)) / 2197 CALL FSUB(N , TK+12*H/13 , UTEMP , V4) DO 25 I=1,N 25 UTEMP(I) = U(I) + H*(8341*V1(I)-32832*V2(I)+29440*V3(I) & -845*V4(I)) / 4104 CALL FSUB(N , TK+H , UTEMP , V5) DO 30 I=1,N 30 UTEMP(I) = U(I) + H*(-6080*V1(I)+41040*V2(I)-28352*V3(I) & +9295*V4(I)-5643*V5(I)) / 20520 CALL FSUB(N , TK+H/2 , UTEMP , V6) E = 0.0 DO 35 I=1,N U4(I) = U(I) + H*(2375*V1(I)+11264*V3(I)+10985*V4(I) & -4104*V5(I)) / 20520 U5(I) = U(I) + H*(33440*V1(I)+146432*V3(I)+142805*V4(I) & -50787*V5(I)+10260*V6(I)) / 282150 35 E = MAX(E , ABS(U4(I)-U5(I))/H) C IF ERROR TOO LARGE, REPEAT STEP WITH H/2 IF (E .GT. TOL) THEN H = H/2 GO TO 5 ENDIF C OTHERWISE, UPDATE TK AND UK TK = TKP1 DO 40 I=1,N U(I) = U5(I) 40 CONTINUE PRINT 45, TK,(U(I),I=1,N) 45 FORMAT (1X,5E15.5,/,(16X,4E15.5)) C IF ERROR MUCH SMALLER THAN REQUIRED, C DOUBLE VALUE OF H FOR NEXT STEP. IF (E .LT. 0.01*TOL) H = 2*H IF (TK.LT.TFINAL) GO TO 5 RETURN END