SUBROUTINE STEP( X, Y, F, NEQN, H, EPS, WT, START, HOLD, K, 1 KOLD, CRASH, PHI, P, YP, PSI, ERRMSG) * * Subroutine STEP integrates a system of first order ordinary differential * equations one step, normally from X to X+H, using a modified divided * difference form of the Adams Pece formulas. Local extrapolation is used to * improve absolute stability and accuracy. The code adjusts its order and step * size to control the local error per unit step in a generalized sense. * Special devices are included to control roundoff error and tyo detect when * the user is requesting too much accuracy. * * More to come!!! * IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL START, CRASH, PHASE1, NORND CHARACTER ERRMSG*78 DIMENSION Y(NEQN), WT(NEQN), PHI(NEQN, 16), P(NEQN), YP(NEQN) DIMENSION PSI(12), ALPHA(12), BETA(12), SIG(13), W(12), V(12) DIMENSION G(13), GSTR(13), TWO(13) EXTERNAL F * ******************************************************************************* * The only machine dependent constants are based on the machine unit roundoff * * error U, which is the smallest positive number such that 1.0+U .GT. 1. The * * user must calculate U with the program MACHIN and insert TWOU = 2.0*U and * * FOURU = 4.0*U in the DATA statment below before calling the code. * ******************************************************************************* * DATA TWOU, FOURU / 2.4D-7, 4.8D-7 / DATA TWO / 2.0D0, 4.0D0, 8.0D0, 1.6D1, 3.2D1, 6.4D1, 1.28D2, 1 2.56D2, 5.12D2, 1.024D3, 2.048D3, 4.096D3, 8.192D3 / DATA GSTR / 5.0D-1, 8.33D-2, 4.17D-2, 2.64D-2, 1.88D-2, 1.43D-2, 1 1.14D-2, 9.36D-3, 7.89D-3, 6.79D-3, 5.92D-3, 5.24D-3, 4.68D-3 / DATA G(1), G(2) / 1.0D0, 0.5D0 /, SIG(1) /1.0D0 / * * **** Begin Block 0 **** * * Check if step size or error tolerance is too small fro machine precision. * If first step, initialize PHI array and estimate a starting step size. * If step too small, determine an acceptable one. * CRASH = .TRUE. IF (DABS(H) .GE. FOURU*DABS(X)) GOTO 5 H = DSIGN(FOURU*DABS(X), H) RETURN 5 P5EPS = 0.5D0 * EPS * * If error tolerance too small, increase it to an acceptable value. * ROUND = 0.0D0 DO 10 L = 1, NEQN 10 ROUND = ROUND + (Y(L)/WT(L))**2 ROUND = TWOU * DSQRT(ROUND) IF (P5EPS .GE. ROUND) GOTO 15 EPS = 2.0D0 * ROUND * (1.0D0 + FOURU) RETURN 15 CRASH = .FALSE. IF (.NOT. START) GOTO 99 * * Initialize and compute appropriate step size for first step. * ERRMSG = 'subroutine STEP, the call after label 15.' CALL F(X, Y, YP, ERRMSG) IF (ERRMSG(1:1) .NE. ' ') RETURN SUM = 0.0D0 DO 20 L = 1, NEQN PHI(L, 1) = YP(L) PHI(L, 2) = 0.0D0 SUM = SUM + (YP(L)/WT(L))**2 20 CONTINUE SUM = DSQRT(SUM) ABSH = DABS(H) IF (EPS .LT. 16.0D0*SUM*H*H) ANSH = 0.25D0*DSQRT(EPS/SUM) H = DSIGN(DMAX1(ABSH, FOURU*DABS(X)), H) HOLD = 0.0D0 K = 1 KOLD = 0 START = .FALSE. PHASE1 = .TRUE. NORND = .TRUE. IF (P5EPS .GT. 1.0D2*ROUND) GOTO 99 NORND = .FALSE. DO 25 L = 1, NEQN 25 PHI(L, 15) = 0.0D0 99 IFAIL = 0 * * **** End of Block 0 **** * * **** Begin Block 1 **** * * Compute coefficients of formulas for this step. Avoid computing those * quantities not changed when step is not changed. * 100 KP1 = K + 1 KP2 = K + 2 KM1 = K - 1 KM2 = K - 2 * * NS is the number of steps taken with size H, including the current one. * When K .LT. NS, no coefficients are changed. * IF (H .NE. HOLD) NS = 0 NS = MIN0(NS+1, KOLD+1) NSP1 = NS + 1 IF (K .LT. NS) GOTO 199 * * Compute those components of ALPHA(*), BETA(*), SIG(*) which are changed. * BETA(NS) = 1.0D0 REALNS = DFLOAT(NS) ALPHA(NS) = 1.0D0/REALNS TEMP1 = H * REALNS SIG(NSP1) = 1.0D0 IF (K .LT. NSP1) GOTO 110 DO 105 I = NSP1, K IM1 = I - 1 TEMP2 = PSI(IM1) PSI(IM1) = TEMP1 BETA(I) = BETA(IM1) * PSI(IM1) / TEMP2 TEMP1 = TEMP2 + H ALPHA(I) = H / TEMP1 REALI = DFLOAT(I) SIG(I+1) = REALI * ALPHA(I) * SIG(I) 105 CONTINUE 110 PSI(K) = TEMP1 * * Compute coefficients G(*). * Initialize V(*) and set W(*). G(2) is set in the data statement above. * IF (NS .GT. 1) GOTO 120 DO 115 IQ = 1, K TEMP3 = DFLOAT(IQ * (IQ + 1)) V(IQ) = 1.0D0 / TEMP3 W(IQ) = V(IQ) 115 CONTINUE GOTO 140 * * If order was raised, update diagonal part of V(*). * 120 IF (K .LE. KOLD) GOTO 130 TEMP4 = DFLOAT(K * KP1) V(K) = 1.0D0 / TEMP4 NSM2 = NS - 2 IF (NSM2 .LT. 1) GOTO 130 DO 125 J = 1, NSM2 I = K - J V(I) = V(I) - ALPHA(J+1)*V(I+1) 125 CONTINUE * * Update V(*) and set W(*). * 130 LIMIT1 = KP1 - NS TEMP5 = ALPHA(NS) DO 135 IQ = 1, LIMIT1 V(IQ) = V(IQ) - TEMP5*V(IQ+1) W(IQ) = V(IQ) 135 CONTINUE G(NSP1) = W(1) * * Compute the G(*) in the work vestor W(*). * 140 NSP2 = NS + 2 IF (KP1 .LT. NSP2) GOTO 199 DO 150 I = NSP2, KP1 LIMIT2 = KP2 - I TEMP6 = ALPHA(I-1) DO 145 IQ = 1, LIMIT2 145 W(IQ) = W(IQ) - TEMP6*W(IQ+1) G(I) = W(1) 150 CONTINUE 199 CONTINUE * * **** End of Block 1 **** * * **** Begin Block 2 **** * * Predict a solution P(*), evaluate derivatives using predicted solution. * Estimate local error at order K and errors at orders K, K-1, and K-2 as if * constant step size were used. * * Change PHI to PHI star. * IF (K .LT. NSP1) GOTO 215 DO 210 I = NSP1, K TEMP1 = BETA(I) DO 205 L = 1, NEQN 205 PHI(L, I) = TEMP1 * PHI(L, I) 210 CONTINUE * * Predict solutions and differences. * 215 DO 220 L = 1, NEQN PHI(L, KP2) = PHI(L, KP1) PHI(L, KP1) = 0.D0 P(L) = 0.D0 220 CONTINUE DO 230 J = 1, K I = KP1 - J IP1 = I + 1 TEMP2 = G(I) DO 225 L = 1, NEQN P(L) = P(L) + TEMP2*PHI(L, I) PHI(L, I) = PHI(L, I) + PHI(L, IP1) 225 CONTINUE 230 CONTINUE IF (NORND) GOTO 240 DO 235 L = 1, NEQN TAU = H * P(L) - PHI(L, 15) P(L) = Y(L) + TAU PHI(L, 16) = (P(L) - Y(L)) - TAU 235 CONTINUE GOTO 250 240 DO 245 L = 1, NEQN 245 P(L) = Y(L) + H*P(L) 250 XOLD = X X = X + H ABSH = DABS(H) ERRMSG = 'subroutine STEP, the call after label 250.' CALL F(X, P, YP, ERRMSG) IF (ERRMSG(1:1) .NE. ' ') THEN DO 251 L = 1, NEQN 251 Y(L) = P(L) RETURN ENDIF * * Estimate errors at orders K, K-1, and K-2. * ERKM2 = 0.0D0 ERKM1 = 0.0D0 ERK = 0.0D0 DO 265 L = 1, NEQN TEMP3 = 1.0D0 / WT(L) TEMP4 = YP(L) - PHI(L, 1) IF (KM2) 265, 260, 255 255 ERKM2 = ERKM2 + ((PHI(L, KM1) + TEMP4)*TEMP3)**2 260 ERKM1 = ERKM1 + ((PHI(L, K) + TEMP4)*TEMP3)**2 265 ERK = ERK + (TEMP4 * TEMP3)**2 IF (KM2) 280, 275, 270 270 ERKM2 = ABSH * SIG(KM1) * GSTR(KM2) * DSQRT(ERKM2) 275 ERKM1 = ABSH * SIG(K) * GSTR(KM1) * DSQRT(ERKM1) 280 TEMP5 = ABSH * DSQRT(ERK) ERR = TEMP5 * (G(K) - G(KP1)) ERK = TEMP5 * SIG(KP1) * GSTR(K) KNEW = K * * Test if order should be lowered. * IF (KM2) 299, 290, 285 285 IF (DMAX1(ERKM1, ERKM2) .LE. ERK) KNEW = KM1 GOTO 299 290 IF (ERKM1 .LE. 0.5D0*ERK) KNEW = KM1 * * Test if step is successful. * 299 IF (ERR .LE. EPS) GOTO 400 * * **** End of Block 2 **** * * **** Begin Block 3 **** * * The step is unsuccessful. Restore X, PHI(*,*), and PSI(*). If third * consecutive failure, set order to one. If step fails more than 3 times, * consider an optimal step size. Double error tolerance and return if * estimated step size is too small for machine precision. * * Restore X, PHI(*,*), and PSI(*). * PHASE1 = .FALSE. X = XOLD DO 310 I = 1, K TEMP1 = 1.0D0 / BETA(I) IP1 = I + 1 DO 305 L = 1, NEQN 305 PHI(L, I) = TEMP1 * (PHI(L, I) - PHI(L, IP1)) 310 CONTINUE IF (K .LT. 2) GOTO 320 DO 315 I = 2, K 315 PSI(I-1) = PSI(I) - H * * On third failure, set order to one. Thereafter, use optimal step size. * 320 IFAIL = IFAIL + 1 TEMP2 = 0.5D0 IF (IFAIL - 3) 335, 330, 325 325 IF (P5EPS .LT. 0.25D0*ERK) TEMP2 = DSQRT(P5EPS/ERK) 330 KNEW = 1 335 H = TEMP2 * H K = KNEW IF (DABS(H) .GE. FOURU*DABS(X)) GOTO 340 CRASH = .TRUE. H = DSIGN(FOURU*DABS(X), H) EPS = EPS + EPS RETURN 340 GOTO 100 * **** End of Block 3 **** * * **** Begin Block 4 **** * * The step is successful. Correct the predicted solution, evaluate the * derivatives using the corrected solution and update the differences. * Determine the best order and step size for next step. * 400 KOLD = K HOLD = H * * Correct and evaluate. * TEMP1 = H * G(KP1) IF (NORND) GOTO 410 DO 405 L = 1, NEQN RHO = TEMP1 * (YP(L) - PHI(L, 1)) - PHI(L, 16) Y(L) = P(L) + RHO PHI(L, 15) = (Y(L) - P(L)) - RHO 405 CONTINUE GOTO 420 410 DO 415 L = 1, NEQN 415 Y(L) = P(L) + TEMP1*(YP(L) - PHI(L, 1)) 420 ERRMSG = 'subroutine STEP, the call at 420 label.' CALL F(X, Y, YP, ERRMSG) IF (ERRMSG(1:1) .NE. ' ') RETURN * * Update differences for next step. * DO 425 L = 1, NEQN PHI(L, KP1) = YP(L) - PHI(L, 1) PHI(L, KP2) = PHI(L, KP1) - PHI(L, KP2) 425 CONTINUE DO 435 I = 1, K DO 430 L = 1, NEQN 430 PHI(L, I) = PHI(L, I) + PHI(L, KP1) 435 CONTINUE * * Estimate error at order K+1 unless: * In first phase when always raise order; * Already decided to lower order; * Step size not constant so estimate unreliable. * ERKP1 = 0.0D0 IF ((KNEW .EQ. KM1) .OR. (K .EQ. 12)) PHASE1 = .FALSE. IF (PHASE1) GOTO 450 IF (KNEW .EQ. KM1) GOTO 455 IF (KP1 .GT. NS) GOTO 460 DO 440 L = 1, NEQN 440 ERKP1 = ERKP1 + (PHI(L, KP2) / WT(L))**2 ERKP1 = ABSH * GSTR(KP1) * DSQRT(ERKP1) * * Using estimated error at order K+1, determine appropriate order for next * step. * IF (K .GT. 1) GOTO 445 IF (ERKP1 .GE. 0.5D0*ERK) GOTO 460 GOTO 450 445 IF (ERKM1 .LE. DMIN1(ERK, ERKP1)) GOTO 455 IF ((ERKP1 .GE. ERK) .OR. (K .EQ. 12)) GOTO 460 * * Here ERKP1 .LT. ERK .LT. DMAX1(ERKM1, ERKM2) else order would have been * lowered in Block 2. Thus order is to be raised. * Raise order. * 450 K = KP1 ERK = ERKP1 GOTO 460 * * Lower order. * 455 K = KM1 ERK = ERKM1 * * With new order determine appropriate step size for next step. * 460 HNEW = H + H IF (PHASE1) GOTO 465 IF (P5EPS .GE. ERK*TWO(K+1)) GOTO 465 HNEW = H IF (P5EPS .GE. ERK) GOTO 465 TEMP2 = DFLOAT(K+1) R = (P5EPS/ERK)**(1.0D0/TEMP2) HNEW = ABSH * DMAX1(0.5D0, DMIN1(0.9D0, R)) HNEW = DSIGN(DMAX1(HNEW, FOURU*DABS(X)),H) 465 H = HNEW RETURN * * **** End of Block 4 **** * END