1 SUBROUTINE DDRIV2 (N,T,Y,F,TOUT,MSTATE,NROOT,EPS,EWT,MINT,WORK, 2 8 LENW,IWORK,LENIW,G) 3C***BEGIN PROLOGUE DDRIV2 4C***DATE WRITTEN 790601 (YYMMDD) 5C***REVISION DATE 850924 (YYMMDD) 6C***CATEGORY NO. I1A2,I1A1B 7C***KEYWORDS ODE,STIFF,ORDINARY DIFFERENTIAL EQUATIONS, 8C INITIAL VALUE PROBLEMS,GEAR'S METHOD, 9C DOUBLE PRECISION 10C***AUTHOR KAHANER, D. K., NATIONAL BUREAU OF STANDARDS, 11C SUTHERLAND, C. D., LOS ALAMOS NATIONAL LABORATORY 12C***PURPOSE The function of DDRIV2 is to solve N ordinary differential 13C equations of the form dY(I)/dT = F(Y(I),T), given the 14C initial conditions Y(I) = YI. The program has options to 15C allow the solution of both stiff and non-stiff differential 16C equations. DDRIV2 uses double precision arithmetic. 17C***DESCRIPTION 18C 19C I. ABSTRACT ....................................................... 20C 21C The function of DDRIV2 is to solve N ordinary differential 22C equations of the form dY(I)/dT = F(Y(I),T), given the initial 23C conditions Y(I) = YI. The program has options to allow the 24C solution of both stiff and non-stiff differential equations. 25C DDRIV2 is to be called once for each output point of T. 26C 27C II. PARAMETERS .................................................... 28C 29C (REMEMBER--To run DDRIV2 correctly in double precision, ALL 30C non-integer arguments in the call sequence, including 31C arrays, MUST be declared double precision.) 32C The user should use parameter names in the call sequence of DDRIV2 33C for those quantities whose value may be altered by DDRIV2. The 34C parameters in the call sequence are: 35C 36C N = (Input) The number of differential equations. 37C 38C T = The independent variable. On input for the first call, T 39C is the initial point. On output, T is the point at which 40C the solution is given. 41C 42C Y = The vector of dependent variables. Y is used as input on 43C the first call, to set the initial values. On output, Y 44C is the computed solution vector. This array Y is passed 45C in the call sequence of the user-provided routines F and 46C G. 47C 48C F = A subroutine supplied by the user. The name must be 49C declared EXTERNAL in the user's calling program. This 50C subroutine is of the form: 51C SUBROUTINE F (N, T, Y, YDOT) 52C DOUBLE PRECISION Y(*), YDOT(*) 53C . 54C . 55C YDOT(1) = ... 56C . 57C . 58C YDOT(N) = ... 59C END (Sample) 60C This computes YDOT = F(Y,T), the right hand side of the 61C differential equations. Here Y is a vector of length at 62C least N. The actual length of Y is determined by the 63C user's declaration in the program which calls DDRIV2. 64C Thus the dimensioning of Y in F, while required by FORTRAN 65C convention, does not actually allocate any storage. When 66C this subroutine is called, the first N components of Y are 67C intermediate approximations to the solution components. 68C The user should not alter these values. Here YDOT is a 69C vector of length N. The user should only compute YDOT(I) 70C for I from 1 to N. 71C 72C TOUT = (Input) The point at which the solution is desired. 73C 74C MSTATE = An integer describing the status of integration. The user 75C must initialize MSTATE to +1 or -1. If MSTATE is 76C positive, the routine will integrate past TOUT and 77C interpolate the solution. This is the most efficient 78C mode. If MSTATE is negative, the routine will adjust its 79C internal step to reach TOUT exactly (useful if a 80C singularity exists beyond TOUT.) The meaning of the 81C magnitude of MSTATE: 82C 1 (Input) Means the first call to the routine. This 83C value must be set by the user. On all subsequent 84C calls the value of MSTATE should be tested by the 85C user. Unless DDRIV2 is to be reinitialized, only the 86C sign of MSTATE may be changed by the user. (As a 87C convenience to the user who may wish to put out the 88C initial conditions, DDRIV2 can be called with 89C MSTATE=+1(-1), and TOUT=T. In this case the program 90C will return with MSTATE unchanged, i.e., 91C MSTATE=+1(-1).) 92C 2 (Output) Means a successful integration. If a normal 93C continuation is desired (i.e., a further integration 94C in the same direction), simply advance TOUT and call 95C again. All other parameters are automatically set. 96C 3 (Output)(Unsuccessful) Means the integrator has taken 97C 1000 steps without reaching TOUT. The user can 98C continue the integration by simply calling DDRIV2 99C again. Other than an error in problem setup, the 100C most likely cause for this condition is trying to 101C integrate a stiff set of equations with the non-stiff 102C integrator option. (See description of MINT below.) 103C 4 (Output)(Unsuccessful) Means too much accuracy has 104C been requested. EPS has been increased to a value 105C the program estimates is appropriate. The user can 106C continue the integration by simply calling DDRIV2 107C again. 108C 5 (Output) A root was found at a point less than TOUT. 109C The user can continue the integration toward TOUT by 110C simply calling DDRIV2 again. 111C 112C NROOT = (Input) The number of equations whose roots are desired. 113C If NROOT is zero, the root search is not active. This 114C option is useful for obtaining output at points which are 115C not known in advance, but depend upon the solution, e.g., 116C when some solution component takes on a specified value. 117C The root search is carried out using the user-written 118C function G (see description of G below.) DDRIV2 attempts 119C to find the value of T at which one of the equations 120C changes sign. DDRIV2 can find at most one root per 121C equation per internal integration step, and will then 122C return the solution either at TOUT or at a root, whichever 123C occurs first in the direction of integration. The index 124C of the equation whose root is being reported is stored in 125C the sixth element of IWORK. 126C NOTE: NROOT is never altered by this program. 127C 128C EPS = On input, the requested relative accuracy in all solution 129C components. EPS = 0 is allowed. On output, the adjusted 130C relative accuracy if the input value was too small. The 131C value of EPS should be set as large as is reasonable, 132C because the amount of work done by DDRIV2 increases as 133C EPS decreases. 134C 135C EWT = (Input) Problem zero, i.e., the smallest physically 136C meaningful value for the solution. This is used inter- 137C nally to compute an array YWT(I) = MAX(ABS(Y(I)), EWT). 138C One step error estimates divided by YWT(I) are kept less 139C than EPS. Setting EWT to zero provides pure relative 140C error control. However, setting EWT smaller than 141C necessary can adversely affect the running time. 142C 143C MINT = (Input) The integration method flag. 144C MINT = 1 Means the Adams methods, and is used for 145C non-stiff problems. 146C MINT = 2 Means the stiff methods of Gear (i.e., the 147C backward differentiation formulas), and is 148C used for stiff problems. 149C MINT = 3 Means the program dynamically selects the 150C Adams methods when the problem is non-stiff 151C and the Gear methods when the problem is 152C stiff. 153C MINT may not be changed without restarting, i.e., setting 154C the magnitude of MSTATE to 1. 155C 156C WORK 157C LENW = (Input) 158C WORK is an array of LENW double precision words used 159C internally for temporary storage. The user must allocate 160C space for this array in the calling program by a statement 161C such as 162C DOUBLE PRECISION WORK(...) 163C The length of WORK should be at least 164C 16*N + 2*NROOT + 204 if MINT is 1, or 165C N*N + 9*N + 2*NROOT + 204 if MINT is 2, or 166C N*N + 16*N + 2*NROOT + 204 if MINT is 3, 167C and LENW should be set to the value used. The contents of 168C WORK should not be disturbed between calls to DDRIV2. 169C 170C IWORK 171C LENIW = (Input) 172C IWORK is an integer array of length LENIW used internally 173C for temporary storage. The user must allocate space for 174C this array in the calling program by a statement such as 175C INTEGER IWORK(...) 176C The length of IWORK should be at least 177C 21 if MINT is 1, or 178C N+21 if MINT is 2 or 3, 179C and LENIW should be set to the value used. The contents 180C of IWORK should not be disturbed between calls to DDRIV2. 181C 182C G = A double precision FORTRAN function supplied by the user 183C if NROOT is not 0. In this case, the name must be 184C declared EXTERNAL in the user's calling program. G is 185C repeatedly called with different values of IROOT to 186C obtain the value of each of the NROOT equations for which 187C a root is desired. G is of the form: 188C DOUBLE PRECISION FUNCTION G (N, T, Y, IROOT) 189C DOUBLE PRECISION Y(*) 190C GO TO (10, ...), IROOT 191C 10 G = ... 192C . 193C . 194C END (Sample) 195C Here, Y is a vector of length at least N, whose first N 196C components are the solution components at the point T. 197C The user should not alter these values. The actual length 198C of Y is determined by the user's declaration in the 199C program which calls DDRIV2. Thus the dimensioning of Y in 200C G, while required by FORTRAN convention, does not actually 201C allocate any storage. 202C 203C***LONG DESCRIPTION 204C 205C III. OTHER COMMUNICATION TO THE USER .............................. 206C 207C A. The solver communicates to the user through the parameters 208C above. In addition it writes diagnostic messages through the 209C standard error handling program XERROR. That program will 210C terminate the user's run if it detects a probable problem setup 211C error, e.g., insufficient storage allocated by the user for the 212C WORK array. Messages are written on the standard error message 213C file. At installations which have this error handling package 214C the user should determine the standard error handling file from 215C the local documentation. Otherwise the short but serviceable 216C routine, XERROR, available with this package, can be used. That 217C program writes on logical unit 6 to transmit messages. A 218C complete description of XERROR is given in the Sandia 219C Laboratories report SAND78-1189 by R. E. Jones. 220C 221C B. The first three elements of WORK and the first five elements of 222C IWORK will contain the following statistical data: 223C AVGH The average step size used. 224C HUSED The step size last used (successfully). 225C AVGORD The average order used. 226C IMXERR The index of the element of the solution vector that 227C contributed most to the last error test. 228C NQUSED The order last used (successfully). 229C NSTEP The number of steps taken. 230C NFE The number of evaluations of the right hand side. 231C NJE The number of evaluations of the Jacobian matrix. 232C 233C IV. REMARKS ....................................................... 234C 235C A. On any return from DDRIV2 all information necessary to continue 236C the calculation is contained in the call sequence parameters, 237C including the work arrays. Thus it is possible to suspend one 238C problem, integrate another, and then return to the first. 239C 240C B. There are user-written routines which are only required by 241C DDRIV3 when certain parameters are set. Thus a message warning 242C of unsatisfied externals may be issued during the load or link 243C phase. This message should never refer to F. This message can 244C be ignored if it refers to G and NROOT is 0. A reference to any 245C other unsatisfied external can be ignored. 246C 247C C. If this package is to be used in an overlay situation, the user 248C must declare in the primary overlay the variables in the call 249C sequence to DDRIV2. 250C 251C V. USAGE .......................................................... 252C 253C PROGRAM SAMPLE 254C EXTERNAL F 255C DOUBLE PRECISION WORK(...), Y(...) See II. for 256C INTEGER IWORK(...) required dimensions for 257C WORK and IWORK 258C OPEN(FILE='TAPE6', UNIT=6, STATUS='NEW') 259C N = ... Number of equations 260C T = ... Initial point 261C DO 10 I = 1,N 262C 10 Y(I) = ... Set initial conditions 263C TOUT = T 264C MSTATE = 1 265C NROOT = 0 266C EPS = ... 267C EWT = ... 268C MINT = 1 269C LENW = ... 270C LENIW = ... 271C 20 CALL DDRIV2 (N, T, Y, F, TOUT, MSTATE, NROOT, EPS, EWT, 272C 8 MINT, WORK, LENW, IWORK, LENIW, G) 273C IF (MSTATE .GT. 2) STOP 274C WRITE(6, 100) TOUT, (Y(I), I=1,N) 275C TOUT = TOUT + 1. 276C IF (TOUT .LE. 10.) GO TO 20 277C 100 FORMAT(...) 278C END (Sample) 279C 280C***REFERENCES GEAR, C. W., "NUMERICAL INITIAL VALUE PROBLEMS IN 281C ORDINARY DIFFERENTIAL EQUATIONS", PRENTICE-HALL, 1971. 282C***ROUTINES CALLED DDRIV3,D1MACH,XERROR 283C***END PROLOGUE DDRIV2 284 EXTERNAL F, JACOBN, FA, G 285 DOUBLE PRECISION EPS, EWT, EWTCOM(1), G, HMAX, D1MACH, T, TOUT, 286 8 WORK(*), Y(*) 287 INTEGER IWORK(*) 288 CHARACTER MSG*81 289 PARAMETER(IMPL = 0, ML = 0, MU = 0, NDE = 0, MXSTEP = 1000) 290C***FIRST EXECUTABLE STATEMENT DDRIV2 291 IF (MINT .LT. 1 .OR. MINT .GT. 3) THEN 292 WRITE(MSG, '(''DDRIV21FE Illegal input. Improper value for '', 293 8 ''the integration method flag,'', I8)') MINT 294 CALL XERROR(MSG, 81, 21, 2) 295 RETURN 296 END IF 297 IF (MSTATE .GE. 0) THEN 298 NSTATE = MSTATE 299 NTASK = 1 300 ELSE 301 NSTATE = - MSTATE 302 NTASK = 3 303 END IF 304 EWTCOM(1) = EWT 305 IF (EWT .NE. 0.D0) THEN 306 IERROR = 3 307 ELSE 308 IERROR = 2 309 END IF 310 IF (MINT .EQ. 1) THEN 311 MITER = 0 312 MXORD = 12 313 ELSE IF (MINT .EQ. 2) THEN 314 MITER = 2 315 MXORD = 5 316 ELSE IF (MINT .EQ. 3) THEN 317 MITER = 2 318 MXORD = 12 319 END IF 320 HMAX = SQRT(D1MACH(2)) 321 CALL DDRIV3 (N, T, Y, F, NSTATE, TOUT, NTASK, NROOT, EPS, EWTCOM, 322 8 IERROR, MINT, MITER, IMPL, ML, MU, MXORD, HMAX, WORK, 323 8 LENW, IWORK, LENIW, JACOBN, FA, NDE, MXSTEP, G) 324 IF (MSTATE .GE. 0) THEN 325 MSTATE = NSTATE 326 ELSE 327 MSTATE = - NSTATE 328 END IF 329 END 330