SLATEC Routines --- DDRIV1 ---
*DECK DDRIV1
SUBROUTINE DDRIV1 (N, T, Y, F, TOUT, MSTATE, EPS, WORK, LENW,
8 IERFLG)
C***BEGIN PROLOGUE DDRIV1
C***PURPOSE The function of DDRIV1 is to solve N (200 or fewer)
C ordinary differential equations of the form
C dY(I)/dT = F(Y(I),T), given the initial conditions
C Y(I) = YI. DDRIV1 uses double precision arithmetic.
C***LIBRARY SLATEC (SDRIVE)
C***CATEGORY I1A2, I1A1B
C***TYPE DOUBLE PRECISION (SDRIV1-S, DDRIV1-D, CDRIV1-C)
C***KEYWORDS DOUBLE PRECISION, GEAR'S METHOD, INITIAL VALUE PROBLEMS,
C ODE, ORDINARY DIFFERENTIAL EQUATIONS, SDRIVE, STIFF
C***AUTHOR Kahaner, D. K., (NIST)
C National Institute of Standards and Technology
C Gaithersburg, MD 20899
C Sutherland, C. D., (LANL)
C Mail Stop D466
C Los Alamos National Laboratory
C Los Alamos, NM 87545
C***DESCRIPTION
C
C Version 92.1
C
C I. CHOOSING THE CORRECT ROUTINE ...................................
C
C SDRIV
C DDRIV
C CDRIV
C These are the generic names for three packages for solving
C initial value problems for ordinary differential equations.
C SDRIV uses single precision arithmetic. DDRIV uses double
C precision arithmetic. CDRIV allows complex-valued
C differential equations, integrated with respect to a single,
C real, independent variable.
C
C As an aid in selecting the proper program, the following is a
C discussion of the important options or restrictions associated with
C each program:
C
C A. DDRIV1 should be tried first for those routine problems with
C no more than 200 differential equations (DDRIV2 and DDRIV3
C have no such restriction.) Internally this routine has two
C important technical defaults:
C 1. Numerical approximation of the Jacobian matrix of the
C right hand side is used.
C 2. The stiff solver option is used.
C Most users of DDRIV1 should not have to concern themselves
C with these details.
C
C B. DDRIV2 should be considered for those problems for which
C DDRIV1 is inadequate. For example, DDRIV1 may have difficulty
C with problems having zero initial conditions and zero
C derivatives. In this case DDRIV2, with an appropriate value
C of the parameter EWT, should perform more efficiently. DDRIV2
C provides three important additional options:
C 1. The nonstiff equation solver (as well as the stiff
C solver) is available.
C 2. The root-finding option is available.
C 3. The program can dynamically select either the non-stiff
C or the stiff methods.
C Internally this routine also defaults to the numerical
C approximation of the Jacobian matrix of the right hand side.
C
C C. DDRIV3 is the most flexible, and hence the most complex, of
C the programs. Its important additional features include:
C 1. The ability to exploit band structure in the Jacobian
C matrix.
C 2. The ability to solve some implicit differential
C equations, i.e., those having the form:
C A(Y,T)*dY/dT = F(Y,T).
C 3. The option of integrating in the one step mode.
C 4. The option of allowing the user to provide a routine
C which computes the analytic Jacobian matrix of the right
C hand side.
C 5. The option of allowing the user to provide a routine
C which does all the matrix algebra associated with
C corrections to the solution components.
C
C II. PARAMETERS ....................................................
C
C (REMEMBER--To run DDRIV1 correctly in double precision, ALL
C non-integer arguments in the call sequence, including
C arrays, MUST be declared double precision.)
C
C The user should use parameter names in the call sequence of DDRIV1
C for those quantities whose value may be altered by DDRIV1. The
C parameters in the call sequence are:
C
C N = (Input) The number of differential equations, N .LE. 200
C
C T = The independent variable. On input for the first call, T
C is the initial point. On output, T is the point at which
C the solution is given.
C
C Y = The vector of dependent variables. Y is used as input on
C the first call, to set the initial values. On output, Y
C is the computed solution vector. This array Y is passed
C in the call sequence of the user-provided routine F. Thus
C parameters required by F can be stored in this array in
C components N+1 and above. (Note: Changes by the user to
C the first N components of this array will take effect only
C after a restart, i.e., after setting MSTATE to +1(-1).)
C
C F = A subroutine supplied by the user. The name must be
C declared EXTERNAL in the user's calling program. This
C subroutine is of the form:
C SUBROUTINE F (N, T, Y, YDOT)
C DOUBLE PRECISION Y(*), YDOT(*)
C .
C .
C YDOT(1) = ...
C .
C .
C YDOT(N) = ...
C END (Sample)
C This computes YDOT = F(Y,T), the right hand side of the
C differential equations. Here Y is a vector of length at
C least N. The actual length of Y is determined by the
C user's declaration in the program which calls DDRIV1.
C Thus the dimensioning of Y in F, while required by FORTRAN
C convention, does not actually allocate any storage. When
C this subroutine is called, the first N components of Y are
C intermediate approximations to the solution components.
C The user should not alter these values. Here YDOT is a
C vector of length N. The user should only compute YDOT(I)
C for I from 1 to N. Normally a return from F passes
C control back to DDRIV1. However, if the user would like
C to abort the calculation, i.e., return control to the
C program which calls DDRIV1, he should set N to zero.
C DDRIV1 will signal this by returning a value of MSTATE
C equal to +5(-5). Altering the value of N in F has no
C effect on the value of N in the call sequence of DDRIV1.
C
C TOUT = (Input) The point at which the solution is desired.
C
C MSTATE = An integer describing the status of integration. The user
C must initialize MSTATE to +1 or -1. If MSTATE is
C positive, the routine will integrate past TOUT and
C interpolate the solution. This is the most efficient
C mode. If MSTATE is negative, the routine will adjust its
C internal step to reach TOUT exactly (useful if a
C singularity exists beyond TOUT.) The meaning of the
C magnitude of MSTATE:
C 1 (Input) Means the first call to the routine. This
C value must be set by the user. On all subsequent
C calls the value of MSTATE should be tested by the
C user. Unless DDRIV1 is to be reinitialized, only the
C sign of MSTATE may be changed by the user. (As a
C convenience to the user who may wish to put out the
C initial conditions, DDRIV1 can be called with
C MSTATE=+1(-1), and TOUT=T. In this case the program
C will return with MSTATE unchanged, i.e.,
C MSTATE=+1(-1).)
C 2 (Output) Means a successful integration. If a normal
C continuation is desired (i.e., a further integration
C in the same direction), simply advance TOUT and call
C again. All other parameters are automatically set.
C 3 (Output)(Unsuccessful) Means the integrator has taken
C 1000 steps without reaching TOUT. The user can
C continue the integration by simply calling DDRIV1
C again.
C 4 (Output)(Unsuccessful) Means too much accuracy has
C been requested. EPS has been increased to a value
C the program estimates is appropriate. The user can
C continue the integration by simply calling DDRIV1
C again.
C 5 (Output)(Unsuccessful) N has been set to zero in
C SUBROUTINE F.
C 6 (Output)(Successful) For MSTATE negative, T is beyond
C TOUT. The solution was obtained by interpolation.
C The user can continue the integration by simply
C advancing TOUT and calling DDRIV1 again.
C 7 (Output)(Unsuccessful) The solution could not be
C obtained. The value of IERFLG (see description
C below) for a "Recoverable" situation indicates the
C type of difficulty encountered: either an illegal
C value for a parameter or an inability to continue the
C solution. For this condition the user should take
C corrective action and reset MSTATE to +1(-1) before
C calling DDRIV1 again. Otherwise the program will
C terminate the run.
C
C EPS = On input, the requested relative accuracy in all solution
C components. On output, the adjusted relative accuracy if
C the input value was too small. The value of EPS should be
C set as large as is reasonable, because the amount of work
C done by DDRIV1 increases as EPS decreases.
C
C WORK
C LENW = (Input)
C WORK is an array of LENW double precision words used
C internally for temporary storage. The user must allocate
C space for this array in the calling program by a statement
C such as
C DOUBLE PRECISION WORK(...)
C The length of WORK should be at least N*N + 11*N + 300
C and LENW should be set to the value used. The contents of
C WORK should not be disturbed between calls to DDRIV1.
C
C IERFLG = An error flag. The error number associated with a
C diagnostic message (see Section IV-A below) is the same as
C the corresponding value of IERFLG. The meaning of IERFLG:
C 0 The routine completed successfully. (No message is
C issued.)
C 3 (Warning) The number of steps required to reach TOUT
C exceeds 1000 .
C 4 (Warning) The value of EPS is too small.
C 11 (Warning) For MSTATE negative, T is beyond TOUT.
C The solution was obtained by interpolation.
C 15 (Warning) The integration step size is below the
C roundoff level of T. (The program issues this
C message as a warning but does not return control to
C the user.)
C 21 (Recoverable) N is greater than 200 .
C 22 (Recoverable) N is not positive.
C 26 (Recoverable) The magnitude of MSTATE is either 0 or
C greater than 7 .
C 27 (Recoverable) EPS is less than zero.
C 32 (Recoverable) Insufficient storage has been allocated
C for the WORK array.
C 41 (Recoverable) The integration step size has gone
C to zero.
C 42 (Recoverable) The integration step size has been
C reduced about 50 times without advancing the
C solution. The problem setup may not be correct.
C 999 (Fatal) The magnitude of MSTATE is 7 .
C
C III. USAGE ........................................................
C
C PROGRAM SAMPLE
C EXTERNAL F
C DOUBLE PRECISION ALFA, EPS, T, TOUT
C C N is the number of equations
C PARAMETER(ALFA = 1.D0, N = 3, LENW = N*N + 11*N + 300)
C DOUBLE PRECISION WORK(LENW), Y(N+1)
C C Initial point
C T = 0.00001D0
C C Set initial conditions
C Y(1) = 10.D0
C Y(2) = 0.D0
C Y(3) = 10.D0
C C Pass parameter
C Y(4) = ALFA
C TOUT = T
C MSTATE = 1
C EPS = .001D0
C 10 CALL DDRIV1 (N, T, Y, F, TOUT, MSTATE, EPS, WORK, LENW,
C 8 IERFLG)
C IF (MSTATE .GT. 2) STOP
C WRITE(*, '(4E12.3)') TOUT, (Y(I), I=1,3)
C TOUT = 10.D0*TOUT
C IF (TOUT .LT. 50.D0) GO TO 10
C END
C
C SUBROUTINE F (N, T, Y, YDOT)
C DOUBLE PRECISION ALFA, T, Y(*), YDOT(*)
C ALFA = Y(N+1)
C YDOT(1) = 1.D0 + ALFA*(Y(2) - Y(1)) - Y(1)*Y(3)
C YDOT(2) = ALFA*(Y(1) - Y(2)) - Y(2)*Y(3)
C YDOT(3) = 1.D0 - Y(3)*(Y(1) + Y(2))
C END
C
C IV. OTHER COMMUNICATION TO THE USER ...............................
C
C A. The solver communicates to the user through the parameters
C above. In addition it writes diagnostic messages through the
C standard error handling program XERMSG. A complete description
C of XERMSG is given in "Guide to the SLATEC Common Mathematical
C Library" by Kirby W. Fong et al.. At installations which do not
C have this error handling package the short but serviceable
C routine, XERMSG, available with this package, can be used. That
C program uses the file named OUTPUT to transmit messages.
C
C B. The number of evaluations of the right hand side can be found
C in the WORK array in the location determined by:
C LENW - (N + 50) + 4
C
C V. REMARKS ........................................................
C
C For other information, see Section IV of the writeup for DDRIV3.
C
C***REFERENCES C. W. Gear, Numerical Initial Value Problems in
C Ordinary Differential Equations, Prentice-Hall, 1971.
C***ROUTINES CALLED DDRIV3, XERMSG
C***REVISION HISTORY (YYMMDD)
C 790601 DATE WRITTEN
C 900329 Initial submission to SLATEC.
C***END PROLOGUE DDRIV1