debdf debdf_r radaua radaua_r rodasa rodasa_r seulexa seulexa_r

seulexa_r() void seulexa_r ( int  n, int  ifcn, double *  t, double  y[], double  tout, double  tend, double *  rtol, double *  atol, int  itol, int  ijac, int  mljac, int  mujac, int  imas, int  mlmas, int  mumas, int  mode, double  work[], int  lwork, int  iwork[], int  liwork, int *  info, double *  tt, double  yy[], double  yyp[], int  ldyypd, double  yypd[], int *  irev  ) Initial value problem of ordinary differential equations (extrapolation method based on the linearly implicit Euler method) (Reverse communication version) Purpose This routine computes a numerical solution of a stiff (or differential algebraic) system of first order ordinary differential equations of the form M * dy/dt = f(t, y), y = y0 at t = t0 where y is the n vector and the system consists of n equations. M is the mass matrix. t0 and y0 are the given initial values of t and y, respectively. This code is the rewritten version of the extrapolation method based on the linearly implicit Euler method code, SEULEX (Reference (1)). This routine is the reverse communication version of seulexa. Parameters [in] n Number of differential equations. (n >= 1) [in] ifcn Specifies whether f(t, y) depends on t. = 0: f(t, y) is independent of t (autonomous). = 1: f(t, y) depends on t (non-autonomous). [in,out] t Independent variable t. This routine integrates from the initial value of t to tend. Depending on the setting of mode, this routine 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 routine 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] ijac Specify how to compute the Jacobian. = 0: Jacobian is computed by finite difference approximation. Never returns with irev = 30. = 1: Jacobian is computed by irev = 30 by user. [in] mljac The lower bandwidth of Jacobian. (0 <= mljac <= n) If mljac = n, Jacobian is stored as full matrix. If mljac < n, Jacobian is stored in band matrix form. [in] mujac The upper bandwidth of Jacobian. (0 <= mujac <= n) If mljac = n, mujac is ignored. [in] imas Specify whether mass matrix M is identity matrix. = 0: M is the identity matrix. Never returns with irev = 40. = 1: M is supplied by irev = 40 by user. [in] mlmas The lower bandwidth of mass matrix M. (0 <= mlmas <= n) If mlmas = n, M is stored in full matrix form. If mlmas < n, M is stored in band matrix form. [in] mumas The upper bandwidth of mass matrix M. (0 <= mumas <= n) If mlmas = n, mumas is ignored. [in] mode Mode of operation. This routine 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 routine 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 routine 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 = 1.0e-6) work[2]: Maximum step size (default = abs(tend - t)) work[3]: Decides whether the Jacobian should be recomputed. (work[3] < 1) (default = 0.001) Increase this value (e.g. to 0.1) when Jacobian evaluations are costly. For small systems, this value should be smaller. work[4], work[5]: Parameters for step size selection. (default: work[4] = 0.1, work[5] = 4) The new step size is chosen subject to the restriction facmin/work[5] <= jth hnew/hold <= 1/facmin, facmin = work[4]^(1/(j - 1)). work[6]: Order selection parameter for decreasing order. (default = 0.7) work[7]: Order selection parameter for increasing order. (default = 0.9) work[8], work[9]: Safety factors for step control algorithm. (default: work[8] = 0.6, work[9] = 0.93) hnew = h*work[9]*(work[8]*tol/err)^(1/(j - 1)). work[11], work[12], work[13], work[14]: Estimated works for a call to f, fjac, LU decompositions, forward-backward substitutions, respectively. (default: work[11] = 1, work[12] = 5, work[13] = 1, work[14] = 1) [out] work[1]: Last step size. [in] lwork Size of array work[]. (n*(ldjac + 2*km + 10) + nm1*(ldmas + lde) + 4*km + n*(km*(km + 1))/2 + 50) where ldjac = nm1 (if full matrix), mljac + mujac + 1 (if band matrix), lde = nm1 (if full matrix), 2*mljac + mujac + 1 (if band matrix), ldmas = nm1 (if full matrix), mlmas + mumas + 1 (if band matrix), 1 (if imas = 0), nm1 = n - m1 (m1 = iwork[8]) 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[0]: Specifies whether to transform the Jacobian matrix to Hessenberg form. It is advantageous for large systems with full Jacobian. Not effective for banded Jacobian nor for the system with IMAS = 1. = 0: Transform Jacobian. = 1: Do not transform Jacobian. iwork[1]: Maximum number of allowed steps. (default = 100000) iwork[2]: Maximum number of columns in the extrapolation table (km). (km >= 3) (default = 12) (if km < 3, default value is used.) iwork[3]: Switch for the step size sequence. (default = 2) = 1: 1, 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, ... = 2: 2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64, ... = 3: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, ... = 4: 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, ... (For other values, default value is used.) iwork[4]: Parameter λ of dense output. (0 or 1) (default = 0) (For other values, default value is used.) iwork[8], iwork[9]: If the differential system has the special structure that y(i)' = y(i + m2) (i = 1, ..., m1) with m1 a multiple of m2, and the remaining equations do not explicitly depend on y'(m1), ..., y'(n - 1), a substantial gain in computer time can be achived by setting the parameters iwork[8] = m1 (> 0) and iwork[9] = m2 (> 0) (m1 + m2 <= n). (default: iwork[8] = 0, iwork[9] = iwork[8]) iwork[12]: Specifies when counters iwork[13] to iwork[19] 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[14]: Number of Jacobian evaluations (either analytically or numerically) iwork[15]: Number of all computed steps. iwork[16]: Number of accepted steps. iwork[17]: Number of rejected steps. iwork[18]: Number of LU decompositions. iwork[19]: Number of forward-backward substitutions. (those needed for step size selection are not counted) [in] liwork Size of array iwork[]. (abs(liwork) >= 2*n + km + 70, 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 routine 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 routine 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 routine again with new tout. = 3: Returned at t = tout to provide intermediate result in mode = 3. It is possible to call this routine again with new tout. = 11: (Error) Maximum number of steps exceeded. [out] tt irev = 1〜20: Refer to yyp[]. irev = 30: Refer to yypd[][]. [out] yy[] Array yy[lyy] (lyy >= n) irev = 1〜20: Refer to yyp[]. irev = 30: Refer to yypd[][]. [in] yyp[] Array yyp[lyyp] (lyyp >= n) irev = 1〜20: Computed derivatives dy/dt = f(t, y) at given t = tt and y = yy[] should be set in yyp[] in the next call. [in] ldyypd irev = 30, 40: Leading dimendion of two dimensional array yypd[][]. (ldyypd >= max(ljac, lmas), Refer to lwork for ljac and lmas.) [in] yypd[][] Array yypd[lyypd][ldyypd] (lyypd >= n) irev = 30: Computed Jacobian at given t = tt and y = yy[] should be set in yypd[][] in the next call in full matrix form if mljac = n, or in band matrix form if mljac < n. irev = 40: Mass matrix should be set in yypd[][] in the next call in full matrix form if mljac = n, or in band matrix form if mljac < n. [in,out] irev Control variable for reverse communication interface. [in] irev should be initialized to zero in the first call. In succeeding calls, irev should not be altered. [out] If irev is not zero, complete the following tasks and call this routine again without changing irev. = 0: Computation finished. Check return code in info. = 1 to 20: Set the computed derivatives at t = tt and y = yy[] in yyp[]. Do not alter any variables other than yyp[]. = 30: Set the computed Jacobian at t = tt and y = yy[] in yypd[][]. Do not alter any variables other than yypd[][]. = 40: Set the mass matrix in yypd[][]. Do not alter any variables other than yypd[][]. Details The band matrix form is illustrated by the following example, when n = 6, ml = 2 and mu = 1. General matrix form: a11 a12 0 0 0 0 a21 a22 a23 0 0 0 a31 a32 a33 a34 0 0 0 a42 a43 a44 a45 0 0 0 a53 a54 a55 a56 0 0 0 a64 a65 a66 Band matrix form: * a12 a23 a34 a45 a56 a11 a22 a33 a44 a55 a66 a21 a32 a43 a54 a65 * a31 a42 a53 a64 * * Array elements marked * are not used by the routine. Details Example for the differential system with special structure (setting parameters iwork[8] and iwork[9]). m1 = m2 = n/2 for second order system p' = v v' = g(p, v) where p and v are vectors of dimension n/2. When m1 > 0, Jacobian and mass matrix must be stored as folows. - Only the non-zero elements (rows m1 to n-1) of the Jacobian have to be stored. - If mljac = n-m1, the Jacobian is full matrix. dfi/dyj = d f(i+m1) / d y(j) (i = 1 to n-m1, j = 1 to n) - If 0 <= mljac < n-m1, the Jacobian is band matrix (m1 = m2 * mm). Only the mm+1 submatrices are stored. df(i-j+mujac+1)/dy(j+k*m2) = d f(i+m1) / d y(j+k*m2) (i = 1 to mljac+mujac+1, j = 1 to m2, k = 0 to mm) mljac is the maximal lower bandwidth of these mm+1 submatrices. mujac is the maximal upper bandwidth of these mm+1 submatrices. If mljac = n-m1, mujac is ignored. - Only the elements of right lower block of dimension n-m1, which is not identity matrix, of mass matrix M are stored. - If mlmas = n-m1, this block is full matrix. am[j-1][i-1] = M(i+m1, j+m1) (i = 1 to n-m1, j = 1 to n-m1). - If 0 <= mlmas < n-m1, this block is band matrix. am[j-1][i-j+mumas] = M(i+m1, j+m1) (i = 1 to n-m1, j = 1 to mlmas+mumas+1). mlmas is the lower bandwidth of the submatrix. mumas is the upper bandwidth of the submatrix. If mlmas = n-m1, mumas is ignored. Reference (1) E. Hairer, S.P. Norsett and G. Wanner, "Solving Ordinary Differential Equations II. Stiff and differential-algebraic Problems. 2nd edition", Springer Series in Computational Mathematics, Springer-Verlag (1996)