DVODERun Method

Namespace: DotNumerics.ODE.DVode
Assembly: DWSIM.MathOps.DotNumerics (in DWSIM.MathOps.DotNumerics.dll) Version: 1.0.0.0 (1.0.0.0)

Syntax

Copy

public void Run(
	IFEX F,
	int NEQ,
	ref double[] Y,
	int offset_y,
	ref double T,
	double TOUT,
	int ITOL,
	double[] RTOL,
	int offset_rtol,
	double[] ATOL,
	int offset_atol,
	int ITASK,
	ref int ISTATE,
	int IOPT,
	ref double[] RWORK,
	int offset_rwork,
	int LRW,
	ref int[] IWORK,
	int offset_iwork,
	int LIW,
	IJEX JAC,
	int MF,
	double[] RPAR,
	int offset_rpar,
	int[] IPAR,
	int offset_ipar
)

Public Sub Run ( 
	F As IFEX,
	NEQ As Integer,
	ByRef Y As Double(),
	offset_y As Integer,
	ByRef T As Double,
	TOUT As Double,
	ITOL As Integer,
	RTOL As Double(),
	offset_rtol As Integer,
	ATOL As Double(),
	offset_atol As Integer,
	ITASK As Integer,
	ByRef ISTATE As Integer,
	IOPT As Integer,
	ByRef RWORK As Double(),
	offset_rwork As Integer,
	LRW As Integer,
	ByRef IWORK As Integer(),
	offset_iwork As Integer,
	LIW As Integer,
	JAC As IJEX,
	MF As Integer,
	RPAR As Double(),
	offset_rpar As Integer,
	IPAR As Integer(),
	offset_ipar As Integer
)

Request Example View Source

Parameters

F IFEX: = The name of the user-supplied subroutine defining the ODE system. The system must be put in the first-order form dy/dt = f(t,y), where f is a vector-valued function of the scalar t and the vector y. Subroutine F is to compute the function f. It is to have the form SUBROUTINE F (NEQ, T, Y, YDOT, RPAR, IPAR) DOUBLE PRECISION T, Y(NEQ), YDOT(NEQ), RPAR where NEQ, T, and Y are input, and the array YDOT = f(t,y) is output. Y and YDOT are arrays of length NEQ. Subroutine F should not alter Y(1),...,Y(NEQ). F must be declared EXTERNAL in the calling program. Subroutine F may access user-defined real and integer work arrays RPAR and IPAR, which are to be dimensioned in the main program. If quantities computed in the F routine are needed externally to DVODE, an extra call to F should be made for this purpose, for consistent and accurate results. If only the derivative dy/dt is needed, use DVINDY instead.
NEQ Int32: = The size of the ODE system (number of first order ordinary differential equations). Used only for input. NEQ may not be increased during the problem, but can be decreased (with ISTATE = 3 in the input).
Y Double: = A real array for the vector of dependent variables, of length NEQ or more. Used for both input and output on the first call (ISTATE = 1), and only for output on other calls. On the first call, Y must contain the vector of initial values. In the output, Y contains the computed solution evaluated at T. If desired, the Y array may be used for other purposes between calls to the solver. This array is passed as the Y argument in all calls to F and JAC.
offset_y Int32
T Double: = The independent variable. In the input, T is used only on the first call, as the initial point of the integration. In the output, after each call, T is the value at which a computed solution Y is evaluated (usually the same as TOUT). On an error return, T is the farthest point reached.
TOUT Double: = The next value of t at which a computed solution is desired. Used only for input. When starting the problem (ISTATE = 1), TOUT may be equal to T for one call, then should .ne. T for the next call. For the initial T, an input value of TOUT .ne. T is used in order to determine the direction of the integration (i.e. the algebraic sign of the step sizes) and the rough scale of the problem. Integration in either direction (forward or backward in t) is permitted. If ITASK = 2 or 5 (one-step modes), TOUT is ignored after the first call (i.e. the first call with TOUT .ne. T). Otherwise, TOUT is required on every call. If ITASK = 1, 3, or 4, the values of TOUT need not be monotone, but a value of TOUT which backs up is limited to the current internal t interval, whose endpoints are TCUR - HU and TCUR. (See optional output, below, for TCUR and HU.)
ITOL Int32: = An indicator for the type of error control. See description below under ATOL. Used only for input.
RTOL Double: = A relative error tolerance parameter, either a scalar or an array of length NEQ. See description below under ATOL. Input only.
offset_rtol Int32
ATOL Double: = An absolute error tolerance parameter, either a scalar or an array of length NEQ. Input only. The input parameters ITOL, RTOL, and ATOL determine the error control performed by the solver. The solver will control the vector e = (e(i)) of estimated local errors in Y, according to an inequality of the form rms-norm of ( e(i)/EWT(i) ) .le. 1, where EWT(i) = RTOL(i)*abs(Y(i)) + ATOL(i), and the rms-norm (root-mean-square norm) here is rms-norm(v) = sqrt(sum v(i)**2 / NEQ). Here EWT = (EWT(i)) is a vector of weights which must always be positive, and the values of RTOL and ATOL should all be non-negative. The following table gives the types (scalar/array) of RTOL and ATOL, and the corresponding form of EWT(i). ITOL RTOL ATOL EWT(i) 1 scalar scalar RTOL*ABS(Y(i)) + ATOL 2 scalar array RTOL*ABS(Y(i)) + ATOL(i) 3 array scalar RTOL(i)*ABS(Y(i)) + ATOL 4 array array RTOL(i)*ABS(Y(i)) + ATOL(i) When either of these parameters is a scalar, it need not be dimensioned in the user's calling program. If none of the above choices (with ITOL, RTOL, and ATOL fixed throughout the problem) is suitable, more general error controls can be obtained by substituting user-supplied routines for the setting of EWT and/or for the norm calculation. See Part iv below. If global errors are to be estimated by making a repeated run on the same problem with smaller tolerances, then all components of RTOL and ATOL (i.e. of EWT) should be scaled down uniformly.
offset_atol Int32
ITASK Int32: = An index specifying the task to be performed. Input only. ITASK has the following values and meanings. 1 means normal computation of output values of y(t) at t = TOUT (by overshooting and interpolating). 2 means take one step only and return. 3 means stop at the first internal mesh point at or beyond t = TOUT and return. 4 means normal computation of output values of y(t) at t = TOUT but without overshooting t = TCRIT. TCRIT must be input as RWORK(1). TCRIT may be equal to or beyond TOUT, but not behind it in the direction of integration. This option is useful if the problem has a singularity at or beyond t = TCRIT. 5 means take one step, without passing TCRIT, and return. TCRIT must be input as RWORK(1). Note: If ITASK = 4 or 5 and the solver reaches TCRIT (within roundoff), it will return T = TCRIT (exactly) to indicate this (unless ITASK = 4 and TOUT comes before TCRIT, in which case answers at T = TOUT are returned first).
ISTATE Int32: = an index used for input and output to specify the the state of the calculation. In the input, the values of ISTATE are as follows. 1 means this is the first call for the problem (initializations will be done). See note below. 2 means this is not the first call, and the calculation is to continue normally, with no change in any input parameters except possibly TOUT and ITASK. (If ITOL, RTOL, and/or ATOL are changed between calls with ISTATE = 2, the new values will be used but not tested for legality.) 3 means this is not the first call, and the calculation is to continue normally, but with a change in input parameters other than TOUT and ITASK. Changes are allowed in NEQ, ITOL, RTOL, ATOL, IOPT, LRW, LIW, MF, ML, MU, and any of the optional input except H0. (See IWORK description for ML and MU.) Note: A preliminary call with TOUT = T is not counted as a first call here, as no initialization or checking of input is done. (Such a call is sometimes useful to include the initial conditions in the output.) Thus the first call for which TOUT .ne. T requires ISTATE = 1 in the input. In the output, ISTATE has the following values and meanings. 1 means nothing was done, as TOUT was equal to T with ISTATE = 1 in the input. 2 means the integration was performed successfully. -1 means an excessive amount of work (more than MXSTEP steps) was done on this call, before completing the requested task, but the integration was otherwise successful as far as T. (MXSTEP is an optional input and is normally 500.) To continue, the user may simply reset ISTATE to a value .gt. 1 and call again. (The excess work step counter will be reset to 0.) In addition, the user may increase MXSTEP to avoid this error return. (See optional input below.) -2 means too much accuracy was requested for the precision of the machine being used. This was detected before completing the requested task, but the integration was successful as far as T. To continue, the tolerance parameters must be reset, and ISTATE must be set to 3. The optional output TOLSF may be used for this purpose. (Note: If this condition is detected before taking any steps, then an illegal input return (ISTATE = -3) occurs instead.) -3 means illegal input was detected, before taking any integration steps. See written message for details. Note: If the solver detects an infinite loop of calls to the solver with illegal input, it will cause the run to stop. -4 means there were repeated error test failures on one attempted step, before completing the requested task, but the integration was successful as far as T. The problem may have a singularity, or the input may be inappropriate. -5 means there were repeated convergence test failures on one attempted step, before completing the requested task, but the integration was successful as far as T. This may be caused by an inaccurate Jacobian matrix, if one is being used. -6 means EWT(i) became zero for some i during the integration. Pure relative error control (ATOL(i)=0.0) was requested on a variable which has now vanished. The integration was successful as far as T. Note: Since the normal output value of ISTATE is 2, it does not need to be reset for normal continuation. Also, since a negative input value of ISTATE will be regarded as illegal, a negative output value requires the user to change it, and possibly other input, before calling the solver again.
IOPT Int32: = An integer flag to specify whether or not any optional input is being used on this call. Input only. The optional input is listed separately below. IOPT = 0 means no optional input is being used. Default values will be used in all cases. IOPT = 1 means optional input is being used.
RWORK Double: = A real working array (double precision). The length of RWORK must be at least 20 + NYH*(MAXORD + 1) + 3*NEQ + LWM where NYH = the initial value of NEQ, MAXORD = 12 (if METH = 1) or 5 (if METH = 2) (unless a smaller value is given as an optional input), LWM = length of work space for matrix-related data: LWM = 0 if MITER = 0, LWM = 2*NEQ**2 + 2 if MITER = 1 or 2, and MF.gt.0, LWM = NEQ**2 + 2 if MITER = 1 or 2, and MF.lt.0, LWM = NEQ + 2 if MITER = 3, LWM = (3*ML+2*MU+2)*NEQ + 2 if MITER = 4 or 5, and MF.gt.0, LWM = (2*ML+MU+1)*NEQ + 2 if MITER = 4 or 5, and MF.lt.0. (See the MF description for METH and MITER.) Thus if MAXORD has its default value and NEQ is constant, this length is: 20 + 16*NEQ for MF = 10, 22 + 16*NEQ + 2*NEQ**2 for MF = 11 or 12, 22 + 16*NEQ + NEQ**2 for MF = -11 or -12, 22 + 17*NEQ for MF = 13, 22 + 18*NEQ + (3*ML+2*MU)*NEQ for MF = 14 or 15, 22 + 17*NEQ + (2*ML+MU)*NEQ for MF = -14 or -15, 20 + 9*NEQ for MF = 20, 22 + 9*NEQ + 2*NEQ**2 for MF = 21 or 22, 22 + 9*NEQ + NEQ**2 for MF = -21 or -22, 22 + 10*NEQ for MF = 23, 22 + 11*NEQ + (3*ML+2*MU)*NEQ for MF = 24 or 25. 22 + 10*NEQ + (2*ML+MU)*NEQ for MF = -24 or -25. The first 20 words of RWORK are reserved for conditional and optional input and optional output. The following word in RWORK is a conditional input: RWORK(1) = TCRIT = critical value of t which the solver is not to overshoot. Required if ITASK is 4 or 5, and ignored otherwise. (See ITASK.)
offset_rwork Int32
LRW Int32: = The length of the array RWORK, as declared by the user. (This will be checked by the solver.)
IWORK Int32: = An integer work array. The length of IWORK must be at least 30 if MITER = 0 or 3 (MF = 10, 13, 20, 23), or 30 + NEQ otherwise (abs(MF) = 11,12,14,15,21,22,24,25). The first 30 words of IWORK are reserved for conditional and optional input and optional output. The following 2 words in IWORK are conditional input: IWORK(1) = ML These are the lower and upper IWORK(2) = MU half-bandwidths, respectively, of the banded Jacobian, excluding the main diagonal. The band is defined by the matrix locations (i,j) with i-ML .le. j .le. i+MU. ML and MU must satisfy 0 .le. ML,MU .le. NEQ-1. These are required if MITER is 4 or 5, and ignored otherwise. ML and MU may in fact be the band parameters for a matrix to which df/dy is only approximately equal.
offset_iwork Int32
LIW Int32: = the length of the array IWORK, as declared by the user. (This will be checked by the solver.)
JAC IJEX: = The name of the user-supplied routine (MITER = 1 or 4) to compute the Jacobian matrix, df/dy, as a function of the scalar t and the vector y. It is to have the form SUBROUTINE JAC (NEQ, T, Y, ML, MU, PD, NROWPD, RPAR, IPAR) DOUBLE PRECISION T, Y(NEQ), PD(NROWPD,NEQ), RPAR where NEQ, T, Y, ML, MU, and NROWPD are input and the array PD is to be loaded with partial derivatives (elements of the Jacobian matrix) in the output. PD must be given a first dimension of NROWPD. T and Y have the same meaning as in Subroutine F. In the full matrix case (MITER = 1), ML and MU are ignored, and the Jacobian is to be loaded into PD in columnwise manner, with df(i)/dy(j) loaded into PD(i,j). In the band matrix case (MITER = 4), the elements within the band are to be loaded into PD in columnwise manner, with diagonal lines of df/dy loaded into the rows of PD. Thus df(i)/dy(j) is to be loaded into PD(i-j+MU+1,j). ML and MU are the half-bandwidth parameters. (See IWORK). The locations in PD in the two triangular areas which correspond to nonexistent matrix elements can be ignored or loaded arbitrarily, as they are overwritten by DVODE. JAC need not provide df/dy exactly. A crude approximation (possibly with a smaller bandwidth) will do. In either case, PD is preset to zero by the solver, so that only the nonzero elements need be loaded by JAC. Each call to JAC is preceded by a call to F with the same arguments NEQ, T, and Y. Thus to gain some efficiency, intermediate quantities shared by both calculations may be saved in a user COMMON block by F and not recomputed by JAC, if desired. Also, JAC may alter the Y array, if desired. JAC must be declared external in the calling program. Subroutine JAC may access user-defined real and integer work arrays, RPAR and IPAR, whose dimensions are set by the user in the main program.
MF Int32: = The method flag. Used only for input. The legal values of MF are 10, 11, 12, 13, 14, 15, 20, 21, 22, 23, 24, 25, -11, -12, -14, -15, -21, -22, -24, -25. MF is a signed two-digit integer, MF = JSV*(10*METH + MITER). JSV = SIGN(MF) indicates the Jacobian-saving strategy: JSV = 1 means a copy of the Jacobian is saved for reuse in the corrector iteration algorithm. JSV = -1 means a copy of the Jacobian is not saved (valid only for MITER = 1, 2, 4, or 5). METH indicates the basic linear multistep method: METH = 1 means the implicit Adams method. METH = 2 means the method based on backward differentiation formulas (BDF-s). MITER indicates the corrector iteration method: MITER = 0 means functional iteration (no Jacobian matrix is involved). MITER = 1 means chord iteration with a user-supplied full (NEQ by NEQ) Jacobian. MITER = 2 means chord iteration with an internally generated (difference quotient) full Jacobian (using NEQ extra calls to F per df/dy value). MITER = 3 means chord iteration with an internally generated diagonal Jacobian approximation (using 1 extra call to F per df/dy evaluation). MITER = 4 means chord iteration with a user-supplied banded Jacobian. MITER = 5 means chord iteration with an internally generated banded Jacobian (using ML+MU+1 extra calls to F per df/dy evaluation). If MITER = 1 or 4, the user must supply a subroutine JAC (the name is arbitrary) as described above under JAC. For other values of MITER, a dummy argument can be used.
RPAR Double: User-specified array used to communicate real parameters to user-supplied subroutines. If RPAR is a vector, then it must be dimensioned in the user's main program. If it is unused or it is a scalar, then it need not be dimensioned.
offset_rpar Int32
IPAR Int32: User-specified array used to communicate integer parameter to user-supplied subroutines. The comments on dimensioning RPAR apply to IPAR.
offset_ipar Int32

Reference

DVODE Class

DotNumerics.ODE.DVode Namespace