SUBROUTINE DE(F, NEQN, Y, T, TOUT, RELERR, ABSERR, IFLAG, ERRMSG) * * Subroutine DE integrates a system of up to 20 first order ordinary * differential equations of the form * DY(I)/DT = F(T, Y(1), Y(2), ..., Y(NEQN)) * Y(I) given at T. * The subroutine integrates from T to TOUT. On return the parameters in the * call list are initialized for continuing the integration. The user has * only to define a new value TOUT and call DE again. * * More to come!!! * IMPLICIT DOUBLE PRECISION (A-H, O-Z) LOGICAL START, CRASH, STIFF CHARACTER ERRMSG*78 DIMENSION Y(NEQN), PSI(12) DIMENSION YY(20), WT(20), PHI(20,16), P(20), YP(20), YPOUT(20) EXTERNAL F * ***************************************************************************** * The only machine dependent constant is based on the machine unit roundoff * * error U which is the smallest positive number such that 1.0+U .GT. 1.0. * * U must be calculated and FOURU = 4.0*U inserted in the following data * * statement before using DE. The routine MACHIN calculates U. FOURU and * * TWOU=2.0*U must also be inserted in subroutine STEP before calling DE. * * For the Alpha-class workstations running OSF/1, U = 0.00000012. * ***************************************************************************** * DATA FOURU/4.8D-7/ * * The constant MAXNUM is the maximum number of steps allowed in one call to * DE. The user may change this limit by altering the following statement: * DATA MAXNUM/500/ * * Test for improper parameters. * IF ((NEQN .LT. 1) .OR. (NEQN .GT. 20)) GO TO 10 IF (T .EQ. TOUT) GO TO 10 IF ((RELERR .LT. 0.0D0) .OR. (ABSERR .LT. 0.0D0)) GO TO 10 EPS = DMAX1(RELERR, ABSERR) IF (EPS .LE. 0.0D0) GO TO 10 IF (IFLAG .EQ. 0) GO TO 10 ISN = ISIGN(1, IFLAG) IFLAG = IABS(IFLAG) IF (IFLAG .EQ. 1) GO TO 20 IF (T .NE. TOLD) GO TO 10 IF ((IFLAG .GE. 2) .AND. (IFLAG .LE. 5)) GO TO 20 10 IFLAG = 6 RETURN * 20 DEL = TOUT - T ABSDEL = DABS(DEL) TEND = T + 10.0D0*DEL IF (ISN .LT. 0) TEND = TOUT NOSTEP = 0 KLE4 = 0 STIFF = .FALSE. RELEPS = RELERR/EPS ABSEPS = ABSERR/EPS * write(6, 999) releps, abseps, eps * 999 format(' releps, abseps, eps = ', 1p3e13.4) IF (IFLAG .EQ. 1) GO TO 30 IF (ISNOLD .LT. 0) GO TO 30 IF (DELSGN*DEL .GT. 0.0D0) GO TO 50 * * On start and restart also set work variables X and YY(*), store the * direction of integration and initialize the step size. * 30 START = .TRUE. X = T DO 40 L = 1, NEQN 40 YY(L) = Y(L) DELSGN = DSIGN(1.0D0, DEL) H = DSIGN(DMAX1(DABS(TOUT-X), FOURU*DABS(X)), TOUT-X) * * If already past output point, interpolate and return. * 50 IF (DABS(X-T) .LT. ABSDEL) GO TO 60 CALL INTRP(X, YY, TOUT, Y, YPOUT, NEQN, KOLD, PHI, PSI) IFLAG = 2 T = TOUT TOLD = T ISNOLD = ISN RETURN * * If cannot go past output point and sufficiently close, extrapolate and * return. * 60 IF ((ISN .GT. 0) .OR. (DABS(TOUT-X) .GE. FOURU*DABS(X))) GO TO 80 H = TOUT-X * ERRMSG = 'subroutine DE, the call after label 60.' CALL F(X, YY, YP, ERRMSG) DO 70 L = 1, NEQN 70 Y(L) = YY(L) + H*YP(L) IF (ERRMSG(1:1) .NE. ' ') RETURN * IFLAG = 2 T = TOUT TOLD = T ISNOLD = ISN RETURN * * Test for too much work. * 80 IF (NOSTEP .LT. MAXNUM) GO TO 100 IFLAG = ISN*4 IF (STIFF) IFLAG = ISN*5 DO 90 L =1, NEQN 90 Y(L) = YY(L) T = X TOLD = T ISNOLD = 1 RETURN * * Limit step size, set weight vector and take a step. * 100 H = DSIGN(DMIN1(DABS(H), DABS(TEND-X)), H) DO 110 L = 1, NEQN 110 WT(L) = RELEPS*DABS(YY(L)) + ABSEPS CALL STEP(X, YY, F, NEQN, H, EPS, WT, START, HOLD, K, KOLD, 1 CRASH, PHI, P, YP, PSI, ERRMSG) * * We don't check for any error messages here, since any errors in subroutine * F means that in the boundary conditions need to be modified. Hopefully the * driver program will make the appropriate corrections. * * Test for tolerances too small. * IF (.NOT. CRASH) GO TO 130 IFLAG = ISN*3 RELERR = EPS*RELEPS ABSERR = EPS*ABSEPS DO 120 L = 1, NEQN 120 Y(L) = YY(L) T = X TOLD = T ISNOLD = 1 RETURN * * Augment counter on work and test for stiffness. * 130 NOSTEP = NOSTEP + 1 KLE4 = KLE4 + 1 IF (KOLD .GT. 4) KLE4 = 0 IF (KLE4 .GE. 50) STIFF = .TRUE. GO TO 50 END