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