dassl dassl_r radau radau5 radau5_r radau_r radaup radaup_r rodas rodas_r seulex seulex_r

seulex_r() void seulex_r ( int  n, int  ifcn, double *  t, double  y[], double  tout, double *  rtol, double *  atol, int  itol, int  ijac, int  mljac, int  mujac, int  imas, int  mlmas, int  mumas, 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  ldyypd, double  yypd[], int *  irtrn, int *  irev  ) Initial value problem of ordinary differential equations (extrapolation method based on the linearly implicit Euler method) (reverse communication version) NOTE - THIS ROUTINE IS DEPRECATED AND WILL BE REMOVED IN THE NEXT 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 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. M is the mass matrix. The system can be linearly implicit (M is not I (identity matrix)) or explicit (M = I). seulex is the code of extrapolation method based on the linearly implicit Euler method. It is provided with the step control algorithm and the continuous output feature. See for details in the reference below. seulex_r is the reverse communication version of seulex. Parameters [in] n Number of differential equations. (n >= 1) [in] ifcn Gives information on f. = 0: f(t, y) independent of t (autonomous). = 1: f(t, y) may depend on t (non-autonomous). (For other values, ifcn = 0 is assumed.) [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] ijac Switch for the computation of the Jacobian. = 0: Jacobian is computed internally by finite differences, never return with irev = 30. = 1: To return with irev = 30 to compute Jacobian externally. (For other values, ijac = 0 is assumed.) [in] mljac The lower bandwidth of Jacobian. (0 <= mljac <= n) If mljac = n, Jacobian is stored as n x n full matrix. If mljac < n, Jacobian is stored as band matrix. [in] mujac The upper bandwidth of Jacobian. (0 <= mujac <= n) If mljac = n, mujac is ignored. [in] imas Switch for the mass matrix M. = 0: M is supposed to be the identity matrix, never return with irev = 40. = 1: To return with irev = 40 to supply M externally. (For other values, imas = 0 is assumed.) [in] mlmas The lower bandwidth of mass matrix M. (0 <= mlmas <= n) If mlmas = n, M is stored as n x n full matrix. If mlmas < n, M is stored as band matrix. [in] mumas The upper bandwidth of mass matrix M. (0 <= mumas <= n) If mlmas = n, mumas is ignored. [in] iout Switch to return with irev = 50 or 51 to print out the intermediate solutions. = 0: Never return = 1: Return for intermediate output without dense ouput. = 2: 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[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-3 or 1.0e-5 is good for stiff equations with initial transient. This choice is not very important. (default: = 1.0e-4) work[1]: Maximum step size (default = tout - t) work[2]: Decides whether the Jacobian should be recomputed. Increase to 0.01 for example, when Jacobian evaluations are costly. For small systems it should be smaller. (default = min(0.0001, rtol[0])) 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/(j - 1)). (default: work[3] = 0.1, work[4] = 4) work[5], work[6]: Parameters for the order selection. Step size is decreased if w[k-2] <= w[k-1]*work[5]. Step size is increased if w[k-1] <= w[k-2]*work[6]. (default: work[5] = 0.7, 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.6, work[8] = 0.93) work[10], work[11], work[12], work[13]: Estimated works for a call to fcn, jac, dec, sol, respectively. (default: work[10] = 1, work[11] = 5, work[12] = 1, work[13] = 1) [out] work[0]: Predicted step size of the last accepted step. [in] lwork Size of array work[]. (lwork >= n*(ljac + km + 8) + nm1*(lmas + le) + 4*km + km2*nrdens + 20) where km = iwork[2], km2 = km*(km + 1)/2, nrdens = iwork[7], ljac = nm1 (if mljac = nm1), mljac + mujac + 1 (if mljac < nm1), lmas = 0 (if imas = 0), nm1 (if imas = 1 and mlmas = nm1), mlmas + mumas + 1 (if imas = 1 and mlmas < nm1), le = nm1 (if mljac = nm1), 2*mljac + mujac + 1 (if mljac < nm1), nm1 = n - m1 (m1 = iwork[8]) 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[0]: If iwork[0] != 0, the code transforms the Jacobian matrix to Hessenberg form. This is particularly advantageous for large systems with full Jacobian. It does not work for banded Jacobian (mljac < n) and not for implicit systems (imas = 1). iwork[1]: Maximum number of allowed steps. (default = 100000) (if iwork[1] < 0, default value will be used) iwork[2]: The maximum number of columns in the extrapolation table. (iwork[2] >= 3) (default = 12) (If iwork[2] < 3, default value will be 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 "lambda" of dense output (iwork[4] = 0 or 1) (default = 0) (If iwork[4] != 0 and iwork[4] != 1, 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 0 < nrdens < n, icont[] must be set (see below). iwork[8] and 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). . See below for further details. (default: iwork[8] = 0, iwork[9] = iwork[8]) iwork[12]: Specifies when iwork[13] to iwork[19] are reset to zero. = 0: Reset whenever this routine is called. != 0: Reset if info = 0. [out] iwork[13]: nfcn = number of function evaluations (those for numerical evaluation of the Jacobian are not counted) iwork[14]: njac = number of Jacobian evaluations (either analytically or numerically) 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) iwork[18]: ndec = number of LU decompositions (n-dimensional matrix) iwork[19]: nsol = number of forward-backward substitutions [in] liwork Size of array iwork[]. (liwork >= 2*n + km + 20 (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 >= (km + 2)*nrdens (km: see iwork[2], nrdens: see iwork[7])) [in,out] icont[] Array icont[licont] Component number list for which 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). [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) = -13: The argument mljac or mlmas had an illegal value (mlmas > mljac) = -14: The argument mujac or mumas had an illegal value (mumas > mujac) = -17: The argument lwork had an illegal value (lwork too small) = -18: The argument iwork[8] or iwork[9] had an illegal value (iwork[8] < 0 or iwork[9] < 0 or iwork[8] + iwork[9] > n) = -19: the argument liwork had an illegal value (liwork too small) = -21: The argument lrcont had an illegal value (lrcont too small) = -23: The argument licont had an illegal value (licont too small) = -24: the argument info had an illegal value (info != 0 and info != 1) = -28: The argument ldyypd had an illegal value (ldyypd < max(ljac, lmas)) = -31: The argument irev had an illegal value (irev is invalid) = 1: Successful exit = 2: Interrupted by irtrn (normal return) = 11: Maximum number of steps exceeded [out] tt irev = 1 to 20: 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 20: 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 20: The computed derivatives at given t (= tt) and y (= yy[]), i.e. yyp[i] = dyi/dt = fi(tt, yy[0], ..., yy[n-1]) (i = 0 to n-1), should be given in the next call. [in] ldyypd irev = 30 or 40: Leading dimension of the two dimensional array yypd[][]. (ldyypd >= max(ljac, lmas), see lwork for ljac and lmas) [in] yypd[][] Array yypd[lyypd][ldyypd] (lyypd >= n) irev = 30: Two dimensional Jacobian data should be given in general matrix form if mljac = n, or band matrix form if mljac < n. See below for further details of band matrix form. irev = 40: The mass matrix should be given in general matrix form if mlmas = n, or band matrix form if mlmas < n (when imas = 1). See below for further details of band matrix form. [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 20: User should set the computed derivative values at tt and yy[] in yyp[]. Do not alter any variables other than yyp[]. = 30: User should computes the Jacobian at tt and yy[] analytically and set to yypd[][]. Do not alter any variables other than yypd[][]. = 40: User should provide the mass matrix in yypd[][] (when imas = 1). Do not alter any variables other than yypd[][]. = 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[]. The solution y[i] (0 <= i <= n-1) at the arbitrary point t2 in the interval [tt, t] can be computed by the function call     y[i] = contex(i, t2, rcont, icont); Further 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. Further Details Example for the differential system with special structure (setting parameters iwork[6] and iwork[7]). 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 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)