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

odex() void odex ( int  n, void(*)(int, double, double *, double *)  f, double *  t, double  y[], double  tout, double *  rtol, double *  atol, int  itol, void(*)(int, double, double, double *, int, double *, int *, int *)  solout, int  iout, double  work[], int  lwork, int  iwork[], int  liwork, double  rcont[], int  lrcont, int  icont[], int  licont, int *  info  ) Initial value problem of ordinary differential equations (extrapolation method (GBS algorithm)) NOTE - THIS ROUTINE IS DEPRECATED AND WILL BE REMOVED IN THE NEXT VERSION. Purpose This routine integrates a system of first order ordinary differential equations of the form dy/dt = f(t, y), 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. odex is an extrapolation algorithm (GBS) code, based on the explicit midpoint rule with step size control, order selection and dense output. See for details in the reference below. Parameters [in] n Number of differential equations. (n >= 1) [in] f The user supplied subroutine, which calculates the derivatives of the differential equations, defined as follows. void f(int n, double t, const double y[], double yp[]) { yp[i] = computed derivative at t and y[] (for i = 0 to n-1) } where n is the number of equations, and yp[] is the derivatives at given t and y[], i.e. yp[i] = dyi/dt = fi(t, y[0], ..., y[n-1]) (i = 0 to 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 either forward in t (tout > t) or backward in t (tout < t) is permitted. 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] solout The user supplied subroutine to print out the intermediate solutions, which is called after every successful step, defined as follows. void solout(int nr, double told, double t, double y[], int n, double cont[], int icont[], int *irtrn) { Output the y[] values at nr-th step t. told is the previous value of t. n is the order of equations. The input value of irtrn will be 0, 1 or 2 in the first, intermediate or last call of solout, respectively. irtrn also serves to interrupt the integration. If irtrn is set to the negative value in solout, the integration will be interrupted and exit with info = 2. Dense output is supported by the control information rcont[] and icont[] if iout = 2. The solution y[i] (0 <= i <= n-1) at the arbitrary point t2 in the interval [told, t] can be computed by the function call y[i] = contx1(i, t2, rcont, icont); } [in] iout Switch for calling the subroutine solout. = 0: solout is never called = 1: solout is used for output = 2: solout is used for dense output. In this case, number of components for which dense output is required (nrdens) must be specified in iwork[7]. (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 guess. = 1/(norm of f'), usually 1.0e-1 or 1.0e-3 is good. This choice is not very important. (default: = 1.0e-4) work[1]: Maximum step size (default = tout - t) work[3], work[4]: Parameters for step size selection. The new step size for the j-th diagonal entry is chosen subject to the restriction facmin/work[4] <= j-th hnew/hold <= 1/facmin where facmin = work[3]^(1/(2*j - 1)). (default: work[3] = 0.02, work[4] = 4) work[5], work[6]: Parameters for the order selection. Order is decreased if w[k-2] <= w[k-1]*work[5]. Order is increased if w[k-1] <= w[k-2]*work[6]. (default: work[5] = 0.8, work[6] = 0.9) work[7], work[8]: Safety factors for step control algorithm. hnew = h*work[8]*(work[7]*tol/err)^(1/(j - 1)) (default: work[7] = 0.65, work[8] = 0.94) work[9]: Step size is reduced by this factor if the stability check is negative (default = 0.5) [out] work[0]: Predicted step size of the last accepted step. [in] lwork Size of array work[]. (lwork >= n*(km + 7) + 5*km + 2*km*(km + 1)*nrdens + 20 (km: see iwork[2], nrdens: see iwork[7])) 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 = 10000) (if iwork[1] < 0, default value will be used) iwork[2]: km = maximum number of columns in the extrapolation table. (iwork[2] >= 3) (default = 9) (if iwork[2] < 3, default value will be used) iwork[3]: Switch for the step size sequence. (iwork[3] >= 4 if iout = 2) (if iwork[3] is out of range, default value will be used) = 1: 2, 4, 6, 8, 10, 12, 14, 16, ... = 2: 2, 4, 8, 12, 16, 20, 24, 28, ... = 3: 2, 4, 6, 8, 12, 16, 24, 32, ... = 4: 2, 6, 10, 14, 18, 22, 26, 30, ... = 5: 4, 8, 12, 16, 20, 24, 28, 32, ... (default = 1 if iout = 0 or 1, 4 if iout = 2) iwork[4]: Stability check is avtivated at most iwork[4] times in one line of the extrapolation table (default = 1) iwork[5]: Stability check is avtivated only in the lines 1 to iwork[5] of the extrapolation table (default = 2) iwork[6]: Determines the degree of interpolation formula. mu = 2*kappa - iwork[6] + 1 (1 <= iwork[6] <= 6) (default = 4) (if iwork[6] < 1 or iwork[6] > 6, default value will be used) iwork[7]: nrdens = number of components for which dense output is required. If iwork[7] < 0, nrdens = abs(iwork[7]) is assumed. (default = n) nrdens is effective only if iout = 2. Otherwise, nrdens = 0 is assumed. if nrdens > n, nrdens = n is assumed. If 0 < nrdens < n, icont[] must be set (see below). iwork[8]: Switch for error estimator in the dense output formula. = 0: Activated = 1: Suppressed (default = 0) (effective when iout = 2) 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 >= 2*km + 21 (km: see iwork[2])) 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 dense output. Used only when iout = 2. [in] lrcont Size of array rcont[]. (lrcont >= (2*km + 5)*nrdens (km: see iwork[2], nrdens: see iwork[7])) [in,out] icont[] Array icont[licont] Component number list for which dense output is required. Used only when iout = 2. 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). [out] = -1: The argument n had an illegal value (n < 1) = -6: The argument rtol had an illegal value (rtol < 0 or rtol[i] < 0) = -6: The argument rtol or atol had an illegal value (rtol = 0 and atol = 0, or rtol[i] = 0 and atol[i] = 0) = -7: The argument atol had an illegal value (atol < 0 or atol[i] < 0) = -12: The argument lwork had an illegal value (lwork too small) = -14: The argument liwork had an illegal value (liwork too small) = -16: The argument lrcont had an illegal value (lrcont too small) = -18: The argument licont had an illegal value (licont too small) = -19: the argument info had an illegal value (info != 0 and info != 1) = 1: Successful exit = 2: Interrupted by solout (normal return) = 11: Maximum number of steps exceeded 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)