The VODE solverEquation 3–1 is reliable for the solution of a wide range of stiff initial-value problems. The source code for the solver is provided in the Ansys Chemkin subdirectory source_public\vode.f (PC), source_public/vode.f (UNIX), which is extensively commented as to its implementation.
SUBROUTINE DVODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, 1 ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF, 2 RPAR, IPAR) EXTERNAL F, JAC DOUBLE PRECISION Y, T, TOUT, RTOL, ATOL, RWORK, RPAR INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW, 1 MF, IPAR DIMENSION Y(*), RTOL(*), ATOL(*), RWORK(LRW), IWORK(LIW), 1 RPAR(*), IPAR(*) C----------------------------------------------------------------------- C DVODE.. Variable-coefficient Ordinary Differential Equation solver, C with fixed-leading coefficient implementation. C This version is in double precision. C C DVODE solves the initial value problem for stiff or nonstiff C systems of first order ODEs, C dy/dt = f(t,y) , or, in component form, C dy(i)/dt = f(i) = f(i,t,y(1),y(2),...,y(NEQ)) (i = 1,...,NEQ). C DVODE is a package based on the EPISODE and EPISODEB packages, and C on the ODEPACK user interface standard, with minor modifications. C----------------------------------------------------------------------- C Revision History (YYMMDD) C 890615 Date Written C 890922 Added interrupt/restart ability, minor changes throughout. C 910228 Minor revisions in line format, prologue, etc. C 920227 Modifications by D. Pang: C (1) Applied subgennam to get generic intrinsic names. C (2) Changed intrinsic names to generic in comments. C (3) Added *DECK lines before each routine. C 920721 Names of routines and labeled Common blocks changed, so as C to be unique in combined single/double precision code (ACH). C 920722 Minor revisions to prologue (ACH). C 920831 Conversion to double precision done (ACH). C----------------------------------------------------------------------- C References.. C C 1. P. N. Brown, G. D. Byrne, and A. C. Hindmarsh, "VODE: A Variable C Coefficient ODE Solver," SIAM J. Sci. Stat. Comput., 10 (1989), C pp. 1038-1051. Also, LLNL Report UCRL-98412, June 1988. C 2. G. D. Byrne and A. C. Hindmarsh, "A Polyalgorithm for the C Numerical Solution of Ordinary Differential Equations," C ACM Trans. Math. Software, 1 (1975), pp. 71-96. C 3. A. C. Hindmarsh and G. D. Byrne, "EPISODE: An Effective Package C for the Integration of Systems of Ordinary Differential C Equations," LLNL Report UCID-30112, Rev. 1, April 1977. C 4. G. D. Byrne and A. C. Hindmarsh, "EPISODEB: An Experimental C Package for the Integration of Systems of Ordinary Differential C Equations with Banded Jacobians," LLNL Report UCID-30132, April C 1976. C 5. A. C. Hindmarsh, "ODEPACK, a Systematized Collection of ODE C Solvers," in Scientific Computing, R. S. Stepleman et al., eds., C North-Holland, Amsterdam, 1983, pp. 55-64. C 6. K. R. Jackson and R. Sacks-Davis, "An Alternative Implementation C of Variable Step-Size Multistep Formulas for Stiff ODEs," ACM C Trans. Math. Software, 6 (1980), pp. 295-318. C----------------------------------------------------------------------- C Authors.. C C Peter N. Brown and Alan C. Hindmarsh C Computing and Mathematics Research Division, L-316 C Lawrence Livermore National Laboratory C Livermore, CA 94550 C and C George D. Byrne C Exxon Research and Engineering Co. C Clinton Township C Route 22 East C Annandale, NJ 08801 C----------------------------------------------------------------------- C Summary of usage. C C Communication between the user and the DVODE package, for normal C situations, is summarized here. This summary describes only a subset C of the full set of options available. See the full description for C details, including optional communication, nonstandard options, C and instructions for special situations. See also the example C problem (with program and output) following this summary. C C A. First provide a subroutine of the form.. C C SUBROUTINE F (NEQ, T, Y, YDOT, RPAR, IPAR) C DOUBLE PRECISION T, Y, YDOT, RPAR C DIMENSION Y(NEQ), YDOT(NEQ) C C which supplies the vector function f by loading YDOT(i) with f(i). C C B. Next determine (or guess) whether or not the problem is stiff. C Stiffness occurs when the Jacobian matrix df/dy has an eigenvalue C whose real part is negative and large in magnitude, compared to the C reciprocal of the t span of interest. If the problem is nonstiff, C use a method flag MF = 10. If it is stiff, there are four standard C choices for MF (21, 22, 24, 25), and DVODE requires the Jacobian C matrix in some form. In these cases (MF .gt. 0), DVODE will use a C saved copy of the Jacobian matrix. If this is undesirable because of C storage limitations, set MF to the corresponding negative value C (-21, -22, -24, -25). (See full description of MF below.) C The Jacobian matrix is regarded either as full (MF = 21 or 22), C or banded (MF = 24 or 25). In the banded case, DVODE requires two C half-bandwidth parameters ML and MU. These are, respectively, the C widths of the lower and upper parts of the band, excluding the main C diagonal. Thus the band consists of the locations (i,j) with C i-ML .le. j .le. i+MU, and the full bandwidth is ML+MU+1. C C C. If the problem is stiff, you are encouraged to supply the Jacobian C directly (MF = 21 or 24), but if this is not feasible, DVODE will C compute it internally by difference quotients (MF = 22 or 25). C If you are supplying the Jacobian, provide a subroutine of the form.. C C SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD, RPAR, IPAR) C DOUBLE PRECISION T, Y, PD, RPAR C DIMENSION Y(NEQ), PD(NROWPD,NEQ) C C which supplies df/dy by loading PD as follows.. C For a full Jacobian (MF = 21), load PD(i,j) with df(i)/dy(j), C the partial derivative of f(i) with respect to y(j). (Ignore the C ML and MU arguments in this case.) C For a banded Jacobian (MF = 24), load PD(i-j+MU+1,j) with C df(i)/dy(j), i.e. load the diagonal lines of df/dy into the rows of C PD from the top down. C In either case, only nonzero elements need be loaded. C C D. Write a main program which calls subroutine DVODE once for C each point at which answers are desired. This should also provide C for possible use of logical unit 6 for output of error messages C by DVODE. On the first call to DVODE, supply arguments as follows.. C F = Name of subroutine for right-hand side vector f. C This name must be declared external in calling program. C NEQ = Number of first order ODE-s. C Y = Array of initial values, of length NEQ. C T = The initial value of the independent variable. C TOUT = First point where output is desired (.ne. T). C ITOL = 1 or 2 according as ATOL (below) is a scalar or array. C RTOL = Relative tolerance parameter (scalar). C ATOL = Absolute tolerance parameter (scalar or array). C The estimated local error in Y(i) will be controlled so as C to be roughly less (in magnitude) than C EWT(i) = RTOL*abs(Y(i)) + ATOL if ITOL = 1, or C EWT(i) = RTOL*abs(Y(i)) + ATOL(i) if ITOL = 2. C Thus the local error test passes if, in each component, C either the absolute error is less than ATOL (or ATOL(i)), C or the relative error is less than RTOL. C Use RTOL = 0.0 for pure absolute error control, and C use ATOL = 0.0 (or ATOL(i) = 0.0) for pure relative error C control. Caution.. Actual (global) errors may exceed these C local tolerances, so choose them conservatively. C ITASK = 1 for normal computation of output values of Y at t = TOUT. C ISTATE = Integer flag (input and output). Set ISTATE = 1. C IOPT = 0 to indicate no optional input used. C RWORK = Real work array of length at least.. C 20 + 16*NEQ for MF = 10, C 22 + 9*NEQ + 2*NEQ**2 for MF = 21 or 22, C 22 + 11*NEQ + (3*ML + 2*MU)*NEQ for MF = 24 or 25. C LRW = Declared length of RWORK (in user's DIMENSION statement). C IWORK = Integer work array of length at least.. C 30 for MF = 10, C 30 + NEQ for MF = 21, 22, 24, or 25. C If MF = 24 or 25, input in IWORK(1),IWORK(2) the lower C and upper half-bandwidths ML,MU. C LIW = Declared length of IWORK (in user's DIMENSION). C JAC = Name of subroutine for Jacobian matrix (MF = 21 or 24). C If used, this name must be declared external in calling C program. If not used, pass a placeholder name. C MF = Method flag. Standard values are.. C 10 for nonstiff (Adams) method, no Jacobian used. C 21 for stiff (BDF) method, user-supplied full Jacobian. C 22 for stiff method, internally generated full Jacobian. C 24 for stiff method, user-supplied banded Jacobian. C 25 for stiff method, internally generated banded Jacobian. C RPAR,IPAR = user-defined real and integer arrays passed to F and JAC. C Note that the main program must declare arrays Y, RWORK, IWORK, C and possibly ATOL, RPAR, and IPAR. C C E. The output from the first call (or any call) is.. C Y = Array of computed values of y(t) vector. C T = Corresponding value of independent variable (normally TOUT). C ISTATE = 2 if DVODE was successful, negative otherwise. C -1 means excess work done on this call. (Perhaps wrong MF.) C -2 means excess accuracy requested. (Tolerances too small.) C -3 means illegal input detected. (See printed message.) C -4 means repeated error test failures. (Check all input.) C -5 means repeated convergence failures. (Perhaps bad C Jacobian supplied or wrong choice of MF or tolerances.) C -6 means error weight became zero during problem. (Solution C component i vanished, and ATOL or ATOL(i) = 0.) C F. To continue the integration after a successful return, simply C reset TOUT and call DVODE again. No other parameters need be reset.