deabm deabm_r derkfa derkfa_r dop853a dop853a_r dopn1210 dopn1210_r dopn43 dopn43_r dopn64 dopn64_r dopn86 dopn86_r dopri5a dopri5a_r dverka dverka_r odex2a odex2a_r odexa odexa_r retarda retarda_r ylaga

odexa() void odexa ( int  n, void(*)(int, double, double *, double *)  f, double *  t, double  y[], double  tout, double  tend, double *  rtol, double *  atol, int  itol, int  mode, double  work[], int  lwork, int  iwork[], int  liwork, int *  info  ) Initial value problem of ordinary differential equations (extrapolation method (GBS algorithm)) Purpose This program integrates a system of first order ordinary differential equations of the form dy/dt = f(t, y), y = y0 at t = t0 where y is the n vector and the system consists of n equations. t0 and y0 are the given initial values of t and y, respectively. This program is the rewrite of the extrapolation (GBS algorithm) program ODEX (Reference (1)). 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[]) { Compute dy/dt and store in yp[]. } where n is the number of equations, and dy/dt (= f(t, y[])) is the computed derivative at given t and y[]. [in,out] t Independent variable t. This program integrates from the initial value of t to tend. Depending on the setting of mode, this program may return at tout or every successful step to provide an intermediate result. [in] An initial value of t (t0). [out] t = tend if the integration has been completed. t = tout or t of the latest step depending on mode in the case of intermediate return. [in,out] y[] Array y[ly] (ly >= n) Dependent variable y. [in] Initial value of y (y0) at initial t (t0). [out] The value of y (computed numerical solution) at t. [in] tout Set tout to the point (t) at which an intermediate result is desired in mode = 2 or 3. tout is not referenced in mode = 0 or 1. In mode 2 or 3, y at tout is computed and returned with info = 2 or 3 respectively. To continue the integration to get result at new tout, it is possible to call this program again with new tout and without changing other variables including info. t < tout <= tend is required. However, if the integration is in backward direction, tend <= tout < t is required. If tend is reached before tout, the integration is terminated immediately and t = tend and info = 0 is returned. [in] tend The point at which an integration is completed. If tend is reached, the integration is terminated and t = tend and info = 0 is returned. The integration is possible for both forward (tend > t) and backward (tend < t) direction. [in] rtol Scalar if itol = 0, or array rtol[lrtol] if itol = 1 (lrtol >= n) (rtol or rtol[i] >= 0) The relative error tolerance(s) to specify how accurately compute the solution. This parameter may be a scalar or an array according to the itol. rtol is used together with atol in a local error test at each step. The integration is performed using step sizes which are automatically selected so as to achieve the following equations. (local error) <= rtol*abs(y[i]) + atol (if itol = 0) or (local error) <= rtol[i]*abs(y[i]) + atol[i] (if itol = 1) for each component of y[] (i = 0 to n-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) to specify how accurately compute the solution. This parameter may be a scalar or an array according to the itol. atol is used together with rtol in a local error test at each step (Refer to rtol above). [in] itol Specifies whether the parameters rtol and atol are scalars or arrays. = 0: rtol and atol are scalars (or arrays with one element). = 1: rtol and atol are arrays. [in] mode Mode of operation. This program integrates from t0 to tend. Four modes of operation are provided depending on the timing of returning imtermediate results. When tend is reached, the integration is terminated and info = 0 is returned in any of four modes. = 0: Not return until tend. tout is not referenced. = 1: Returns at every successful step (info = 1 is returned). tout is not referenced. The end point of the last step is returned in t. In the final step, the step size is adjusted to fit into tend. = 2: Returns at tout during the ingegration to provide an intermediate result (t = tout and info = 2 are returned). To continue the integration to get result at new tout, it is possible to call this program again with new tout. In the final step, the step size is adjusted to fit into tout. = 3: Returns at tout during the ingegration to provide an intermediate result (t = tout and info = 3 are returned). To continue the integration to get result at new tout, it is possible to call this program again with new tout. Different from mode = 2, the value of y at tout is computed by interpolation, and the step size is not adjusted even in the last step before tout. The steps through integration are same with those of mode = 1. [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 used. [in] work[0]: Initial step size (default = 0.0001) work[2]: Maximum step size (default = abs(tend - t)) work[4] and work[5]: Parameters for step size selection. The new step size is chosen subject to the restriction facmin/work[5] <= j-th hnew/hold <= 1/facmin where facmin = work[4]^(1/(2*j - 1)). (default work[4] = 0.02, work[5] = 4) work[6] and work[7]: Parameters for the order selection. Order is decreased if w[k-2] <= w[k-1]*work[6]. Order is increased if w[k-1] <= w[k-2]*work[7]. (default: work[6] = 0.8, work[7] = 0.9) work[8] and work[9]: Safety factors for step control algorithm. hnew = h*work[9]*(work[8]*tol/err)^(1/(j-1)). (default: work[8] = 0.65, work[9] = 0.94) work[10]: Step size is reduced by this factor. (default = 0.5) [out] work[1]: Last step size. [in] lwork Size of array work[]. (abs(lwork) >= 12*n + 5*(n + 1)*km + 2*n*km^2 + 40 (km = iwork[2])) If lwork < 0, the abs(lwork) is used, and all parameters work[0] to work[19] are cleared to 0. [in,out] iwork[] Array iwork[liwork] Integer 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 used. [in] iwork[1]: Maximum number of allowed steps. (default = 10000) iwork[2]: Maximum number of columns in the extrapolation table (km). (km >= 3) (default = 9) (If km < 3, default value will be used) iwork[3]: Switch for the step size sequence. (default = 1 (if mode != 3), 4 (if mode = 3)) = 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, ... If mode = 3, iwork[3] must be 4 or 5. If iwork[3] is out of range, default value will be used. 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 (1 <= iwork[6] <= 6). (default = 4) If iwork[6] is out of range, default value will be used. iwork[8]: Switch for error estimator in the dense output formula (effective when mode = 3). (default = 0) = 0: Activated = 1: Suppressed iwork[12]: Specifies when counters iwork[13] to iwork[17] are reset to zero. (default = 0) = 0: Reset if info = 0. = 1: Do not reset even if info = 0. [out] Statistical information (counters). iwork[13]: Number of function evaluations. iwork[15]: Number of all computed steps. iwork[16]: Number of accepted steps. iwork[17]: Number of rejected steps. [in] liwork Size of array iwork[]. (abs(liwork) >= 2*km + 60 (km = iwork[2])) If liwork < 0, the abs(liwork) is used, and all parameters iwork[0] to iwork[19] are cleared to 0. [in,out] info [in] Control code. = 0: Set info = 0 on the initial call to start new problem. All variables will be initialized and the computation will begin. = 1, 2, 3: When returned with info = 1, 2 or 3, it is possible to call this program again with new tout and without changing info to continue the integration. [out] Return code. = 0: Successful exit. Integration to tend completed. < 0: The (-info)-th argument is invalid. = 1: Returned to provide intermediate result in mode = 1. It is possible to call this program again to forward to the next step. = 2: Returned at t = tout to provide intermediate result in mode = 2. It is possible to call this program again with new tout. = 3: Returned at t = tout to provide intermediate result in mode = 3. It is possible to call this program again with new tout. = 11: (Error) Maximum number of steps exceeded. Reference (1) E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations. Nonstiff Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1993)