derkf derkf_int derkf_r dop853 dop853_r dopri5 dopri5_r doprin doprin_r dverk dverk_int dverk_r odex odex2 odex2_r odex_r retard retard_r

retard_r() void retard_r ( int  n, double *  t, double  y[], double  tout, double *  rtol, double *  atol, int  itol, int  iout, double  work[], int  lwork, int  iwork[], int  liwork, double  rcont[], int  lrcont, int  icont[], int  licont, int *  info, double *  tt, double  yy[], double  yyp[], int *  irtrn, int *  irev  ) Initial value problem of delay differential equations (5(4)-th order Dorman-Prince method) (reverse communication version) NOTE - THIS ROUTINE IS DEPRECATED AND WILL BE REMOVED IN THE NEXT VERSION. Purpose This routine integrates a system of first order delay ordinary differential equations of the form dy/dt = f(t, y(t), y(t-τ), ...), y = y0 at t = t0 where t0 and y0 are the given initial values of t and y, respectively. y may be a vector if the above is a system of differential equations. retard is the explicit Runge-Kutta code based on the 5(4)-th order Dormand-Prince method. It is provided with the step control algorithm and the dense output feature. See for details in the reference below. retard_r is the reverse communication version of retard. Parameters [in] n Number of differential equations. (n >= 1) [in,out] t This routine integrates from t to tout. The initial point of the integration is to be given, and the last point of the final step will be returned. [in] Initial value of the independent variable t. [out] Last value of the independent variable t of the final step (normally equals to tout). The solution was successfully advanced to this point. It is possible to continue the integration to new point by recalling this routine with the new tout value. [in,out] y[] Array y[ly] (ly >= n) [in] Initial values of the dependent variables y[] at initial t. [out] Computed solution approximation at last t (normally equals to tout). [in] tout Set tout to the point at which a solution is desired. Integration should be forward in t (tout > t). The routine advances the solution from t to tout using step sizes which are automatically selected so as to achieve the desired accuracy. [in] rtol Scalar if itol = 0, or array rtol[lrtol] if itol = 1 (lrtol >= n) (rtol or all components of rtol[] >= 0) The relative error tolerance(s). The routine keeps, roughly, the local error of y[i] below rtol*abs(y[i]) + atol (i = 0 to n-1) (if itol = 0) or rtol[i]*abs(y[i]) + atol[i] (i = 0 to n-1) (if itol = 1). Both rtol and atol, or, rtol[i] and atol[i] (i = 0 to n-1) should not be 0 at the same time. [in] atol Scalar if itol = 0, or array atol[latol] if itol = 1 (latol >= n) (atol or all components of atol[] >= 0) The absolute error tolerance(s). The routine keeps, roughly, the local error of y[i] below rtol*abs(y[i]) + atol (i = 0 to n-1) (if itol = 0) or rtol[i]*abs(y[i]) + atol[i] (i = 0 to n-1) (if itol = 1). Both rtol and atol, or, rtol[i] and atol[i] (i = 0 to n-1) should not be 0 at the same time. [in] itol Specifies whether the parameters rtol and atol are scalars or arrays. = 0: rtol and atol are scalars. = 1: rtol and atol are arrays. (For other values, itol = 0 is assumed.) [in] iout Switch to return to print out the intermediate solutions (see irev = 50, 51). = 0: Not return = 1: Return for intermediate output with dense output. In this case, number of components for which dense output is required (nrdens) must be specified in iwork[4]. (For other values, iout = 0 is assumed.) [in,out] work[] Array work[lwork] Work array. work[0] to work[19] serve as parameters for the program. If the input parameter is set to 0, the default parameter value defined for each parameter will be loaded. [in] work[0]: Initial step size (default = estimated by the code) work[1]: Maximal step size (default = tout - t) work[3] and work[4]: Parameters for step size selection. The new step size is chosen subject to the restriction work[3] <= hnew/hold <= work[4]. (default work[3] = 0.2, work[4] = 10) work[7]: The safety factor in step size prediction (0.0001 < work[7] < 1). If work[7] <= 0.0001 or work[7] >= 1, default value is used. (default = 0.9) work[8]: The "beta" for the stabilized step size control (work[8] <= 0.2). If work[8] < 0, beta is set to 0. (default = 0.04) work[20], ..., work[20+ngrid-1]: Prescribed points, which the integration method has to take as grid points (t < work[20] < ... < work[20+ngrid-1] = tout). [out] work[0]: Last step size [in] lwork Size of array work[]. (lwork >= 8*n + ngrid + 20) (ngrid: see iwork[5]) If lwork < 0, abs(lwork) will be used and work[0] to work[19] will be initialized to zeros. [in,out] iwork[] Array iwork[liwork] Work array. iwork[0] to iwork[19] serve as parameters for the program. If the input parameter is set to 0, the default parameter value defined for each parameter will be loaded. [in] iwork[1]: Maximum number of allowed steps. (default = 100000) (if iwork[1] < 0, default value will be used) iwork[3]: Test for stiffness is activated after every iwork[3] steps. For negative iwork[3], the stiffness test is not activated. (default = 1000) iwork[5]: ngrid = number of grid points. (default = 1) The solution may have the discontinuities in its derivatives (grid points). In such case, by specifying the points, the precision and the computation time can be improved. If iwork[5] <= 0, ngrid = 1 is assumed. If iwork[5] > 0, work[20], ..., work[20 + ngrid - 1] should be set by user. The last point work[20 + ngrid - 1] must be equal to tout (if iwork[5] = 0, the last point will be automatically set by the code). iwork[7]: nrdens = number of components for which dense output is required (either by f or solout). If iwork[7] < 0, nrdens = abs(iwork[7]) is assumed. (default = n) If 0 < nrdens < n, icont[] must be set (see below). iwork[12]: Specifies when iwork[13] to iwork[17] are reset to zero. = 0: Reset whenever this routine is called. != 0: Reset if info = 0. [out] iwork[13]: nfcn = number of function evaluations iwork[15]: nstep = number of computed steps iwork[16]: naccept = number of accepted steps iwork[17]: nreject = number of rejected steps (due to error test) (step rejections in the first step are not counted) [in] liwork Size of array iwork[]. (liwork >= 20) If liwork < 0, abs(liwork) will be used and iwork[0] to iwork[19] will be initialized to zeros. [out] rcont[] Array rcont[lrcont] Control information area for computing y(t-τ) and dense output. [in] lrcont Size of array rcont[]. (lrcont >= mxst*(5*nrdens + 2) (nrdens: see iwork[7], mxst: maximum number of back steps to be stored (back steps available through ylag))) [in,out] icont[] Array icont[licont] Component number list for which y(t-τ) value and/or dense output is required. If nrdens = n or iwork[7] < 0, it will be automatically set to 1, ..., nrdens. If 0 < nrdens < n, user must set the component numbers for which dense output is required. See iwork[7] for nrdens. [in] licont Size of array icont[]. (licont >= nrdens (nrdens: see iwork[7])) [in,out] info [in] = 0: Initialize and start computation (Solve new problem). = 1: Continue computation with new tout value (Resume computation of previous call). Available only when Ngrod = 1. [out] = -1: The argument n had an illegal value (n < 1) = -5: The argument rtol had an illegal value (rtol < 0 or rtol[i] < 0) = -5: The argument rtol or atol had an illegal value (rtol = 0 and atol = 0, or rtol[i] = 0 and atol[i] = 0) = -6: The argument atol had an illegal value (atol < 0 or atol[i] < 0) = -9: The argument work[20 + i - 1] had an illegal value (not in ascending order) = -10: The argument lwork had an illegal value (lwork too small) = -12: The argument liwork had an illegal value (liwork too small) = -14: The argument lrcont had an illegal value (lrcont too small) = -16: The argument licont had an illegal value (licont too small) = -17: the argument info had an illegal value (info != 0 and info != 1) = -22: the argument irev had an illegal value (irev != 1, ..., 9, 50 nor 51) = 1: Successful exit = 2: Interrupted by irtrn (normal return) = 11: Maximum number of steps exceeded = 12: Step size becomes too small = 13: Problem is probably stiff (interrupted) = 14: Computation interrupted by ylag (e.g. mxst too small) [out] tt irev = 1 to 9: The value of t where the derivative values should be evaluated and given in yyp[] in the next call. [out] yy[] Array yy[lyy] (lyy >= n) irev = 1 to 9: The value of y where the derivative values should be evaluated and given in yyp[] in the next call. [in] yyp[] Array yyp[lyyp] (lyyp >= n) irev = 1 to 9: The computed derivatives at given t (= tt) and y (= yy[]), i.e. yyp[i] = dyi/dt = fi(t, y(t), y(t-τ), ...) (i = 0 to n-1), should be given in the next call. y(t - τ) can be obtained by the reverse communication subroutine call ylag_r(i, x, rcont, icont, &yy, &irev); where i is the component number (0 to n - 1) and x is t - τ. rcont and icont must be passed without altering those values. If irev = 0, the result is returned in yy. If t0 - τ <= x <= t0 (where t0 is the initial value of t), irev = 1 is returned. Then an initial value of y[i] at X must be given in yy. [in,out] irtrn [in] irev = 50 or 51: Do not alter irtrn unless user want to interrupt the integration. If irtrn is set to the negative value, the integration will be interrupted and exit with info = 2. [out] irev = 50 or 51: The output value will be 0, 1 or 2 in the first, intermediate or last return with irev = 50 or 51, respectively. [in,out] irev Control variable for reverse communication. [in] Before first call, irev should be initialized to zero. On succeeding calls, irev should not be altered. [out] If irev is not zero, complete the following process and call this routine again. = 0: Computation finished. See return code in info. = 1 to 9: User should set the computed derivative values at tt and yy[] in yyp[]. Do not alter any variables other than yyp[]. = 50, 51: To be returned with this code to print out the intermediate solutions after every successful step if iout = 1. The values of y[] at t are provided. tt is the previous value of t. Dense output is supported by the control information rcont[] and icont[] if iout = 1. The solution y[i] (0 <= i <= n-1) at the arbitrary point before t can be computed by calling ylag_r(). See description of yyp[]. Reference E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1993)