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

dassl_r() void dassl_r ( int  n, double *  t, double  y[], double  yp[], double  tout, int  iopt[], double *  rtol, double *  atol, double  work[], int  lwork, int  iwork[], int  liwork, int *  info, double *  tt, double  yyp[], int  ldyypd, double  yypd[], double *  cj, int *  ires, int *  irev  ) Solution of differential algebraic equation (DAE) (1~5-th order backward differentiation formula (BDF)) (reverse communication version) NOTE - THIS ROUTINE IS DEPRECATED AND WILL BE REMOVED IN THE NEXT VERSION. Purpose This routine solves a implicit system of differential algebraic equations (DAEs) of the form f(t, y, y') = 0, y = y0 and y' = y'0 at t = t0 where t0, y0 and y'0 are the given initial values of t, y and y', respectively. y and y' may be vectors if the above is a system of differential algebraic equations. In this routine, the derivatives y' are approximated by backward differentiation formulae (BDF) of orders 1 through 5, and the resulting nonlinear system at each time-step is solved by Newton's method. dassl_r is the reverse communication version of dassl. 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. It is possible to continue the integration to get result at new tout. This is the interval mode of operation. t should be equal to the previous tout on continuation call. It is also possible for the routine to return with the solution at each intermediate step on the way to tout. This is the intermediate-output mode of operation. This mode is a good way to proceed if you want to see the behavior of the solution. The mode of operation is specified by the parameter iopt[2]. [in] Initial value of the independent variable t. [out] Last value of the independent variable t of the final step (normally equals to tout in interval mode). The solution was successfully advanced to this point. [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 (= tout in interval mode). [in,out] yp[] Array yp[lyp] (lyp >= n) [in] Initial values of the derivatives y' at initial t. [out] Computed derivatives of the solution approximation at last t (= tout in interval mode). [in] tout Set tout to the point at which a solution is desired. You can not take tout = t. Integration either forward in t (tout > t) or backward in t (tout < t) is permitted. It is, however, not allowed to change the direction of integration without restarting. The code advances the solution from t to tout using step sizes which are automatically selected so as to achieve the desired accuracy. If you wish, the code will return with the solution and its derivative following each intermediate step (intermediate-output mode) so that you can monitor them, but you still must provide tout in accord with the basic aim of the code. [in,out] iopt[] Array iopt[liopt] (liopt >= 15) Use this array to give the code more details about how you want your problem solved. For standard use, set all entries to 0. [in] iopt[0]: This parameter is automatically set by the program. No need to set by user. iopt[1]: Specifies whether the parameters rtol and atol are scalars or arrays. It cannot be changed after the first call for the current problem.   = 0: rtol and atol are scalars.   = 1: rtol and atol are arrays. iopt[2]: Mode of operation.   = 0: Return only at tout (interval mode)   = 1: Return at every step (intermediate output mode) iopt[3]: This routine may integrate past tout and interpolate to obtain the result at tout. However, it may not be permissible to integrate past some specific point because a discontinuity occurs there, or the function or its derivative is not defined beyond that point. In that case, the stopping point tstop can be specified so that the routine will not integrate beyond that point.   = 0: Not define the stopping point   = 1: Define the stopping point by setting work[0] = tstop iopt[4]: If user does not evaluate the partial derivatives of the system of differential equations analytically, it will be approximated by numerical differencing in this code.   = 0: Evaluate the partial derivatives automatically by numerical differences.   = 1: The routine will return with irev = 30 to 33 to evaluate the partial derivatives analytically bu user. iopt[5]: Specify whether the matrix of partial derivatives is computed as full (general) matrix or band matrix. In general, computation with band matrix costs less time and storage than with full matrix if 2*ml + mu < n, where ml is lower bandwidth and mu is upper bandwidth.   = 0: n x n full (general) matrix.   = 1: Band matrix. Bandwidths should be specified by setting iwork[0] = ml and iwork[1] = mu. (0 <= ml <= n and 0 <= mu <= n) iopt[6]: User can specify a maximum (absolute value of) stepsize.   = 0: Maximum step size is decided by the program.   = 1: Specify maximum step size by setting work[1] = hmax. iopt[7]: User can specify an initial (absolute value of) stepsize.   = 0: Initial step size is decided by the program.   = 1: Specify initial step size by setting work[2] = h0. iopt[8]: User can specify the maximum order.   = 0: Default maximum order (maxord = 5) is used.   = 1: Specify maximum order by setting iwork[2] = maxord (1 <= maxord <= 5). iopt[9]: If the solutions to your equations will always be nonnegative, it may help to set this parameter.   = 0: Solve the problem without invoking any special nonnegativity constraints.   = 1: Solve the problem with invoking special nonnegativity constraints. iopt[10]: The initial t, y, and y' are normally consistent, i.e. f(t, y, y') = 0 at the initial time. If you do not know the initial derivative precisely, you can request to try to compute it by the program.   = 0: The initial t, y, and y' are consistent.   = 1: Request to try to compute y' by the program. Set yp[] to the initial approximation values. If you have no idea, set yp[] to 0. [out] iopt[0]: The value of this parameter is necessary for continuation calls. Do not alter. [in,out] 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) to tell the code how accurately you want the solution to be computed. This parameter may be a scalar or an array according to the iopt[1]. The array can be used to control the error test more precisely. The tolerances are used by the code in a local error test at each step which requires roughly that   abs(local error) <= rtol*abs(y[i]) + atol (if itol = 0)     or   abs(local error) <= rtol[i]*abs(y[i]) + atol[i] (if itol = 1) for each component of y[] (i = 0 to n-1). Setting rtol = 0 yields a pure absolute error test on that component. [in] Relative error tolerance(s) for the error test. [out] When returned with info = 4, the error tolarance(s) will be increased to the value appropriate for continuing the computation. Otherwise, it will remain unchanged. [in,out] 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 tell the code how accurately you want the solution to be computed. This parameter may be a scalar or an array according to the iopt[1]. The array can be used to control the error test more precisely. The tolerances are used by the code in a local error test at each step which requires roughly that   abs(local error) <= rtol*abs(y[i]) + atol (if itol = 0)     or   abs(local error) <= rtol[i]*abs(y[i]) + atol[i] (if itol = 1) for each component of y[] (i = 0 to n-1). Setting atol = 0 results in a pure relative error test on that component. [in] Absolute error tolerance(s) for the error test. [out] When returned with info = 4, the error tolarance(s) will be increased to the value appropriate for continuing the computation. Otherwise, it will remain unchanged. [in,out] work[] Array work[lwork] Work array. Some entries serve as parameters for the program. [in] work[0]: If iopt[3] = 1, the stopping point is specified as work[0] = tstop. work[1]: If iopt[6] = 1, the maximum step size is specified as work[1] = hmax. work[2]: If iopt[7] = 1, the initial step size is specified as work[2] = h0. [out] work[2]: The step size to be attempted on the next step. work[3]: The current value of the independent variable. This will be different from t when interpolation has been performed even info = 1. work[6]: The step size used on the last successful step. [in] lwork The length of work[]. lwork >= n^2 + (maxord + 4)*n + 40 (if full matrix of partial derivatives), lwork >= (2*ml + mu + 1)*n + (maxord + 4)*n + 40 + 2*(n/(ml + mu + 1) + 1) (if banded matrix of partial derivatives). See iopt[5] for ml and mu, iopt[8] for maxord. [in,out] iwork[] Array iwork[liwork] Work array. Some entries serve as parameters for the program. [in] iwork[0] and iwork[1]: If iopt[5] = 1, lower and upper bandwidth of the matrix of partial derivatives are specified as iwork[0] = ml and iwork[1] = mu. iwork[2]: If iopt[8] = 1, the maximum order is specified as iwork[2] = maxord. [out] iwork[6]: The order of the method to be attempted on the next step. iwork[7]: The order of the method used on the last step. iwork[10]: The number of steps taken. iwork[11]: The number of calls to res. iwork[12]: The number of evaluations of the matrix of partial derivatives needed. iwork[13]: The total number of error test failures. iwork[14]: The total number of convergence test failures (includes singular iteration matrix failures). [in] liwork The length of iwork[]. (liwork >= n + 20) [out] info [in] Control code. = 0: Set info = 0 on the initial call for the problem (to start new problem). The routine will be initialized and the computation for the new problem will be started. = 1, 2 or 11 to 13: When returned from the routine with info = 1, 2 or 11 to 13, user can reenter without changing info to continue computation. See descriptions below. [out] Return code. By examining this code, user can call this routine again when info = 1, 2 or 11 to 13 as a next action if necessary. = -1: The argument n had an illegal value (n < 1) = -5: The argument tout had an illegal value (tout = t) = -5: The argument tout had an illegal value (tout is behind t) = -5: The argument tout had an illegal value (tout is too close to t) = -6: The argument iopt[] had an illegal value (some element is not zero or one) = -7: The argument rtol and atol had an illegal value (all elements of rtol and atol are zero) = -7: The argument rtol and atol had an illegal value (some element of wt <= 0) = -7: The argument rtol had an illegal value (rtol < 0 or rtol[i] < 0) = -8: The argument atol had an illegal value (atol < 0 or atol[i] < 0) = -9: The argument work[0] had an illegal value (iopt[3] = 1 and tstop is behind t) = -9: The argument work[0] had an illegal value (iopt[3] = 1 and tstop is behind tout) = -9: The argument work[1] had an illegal value (hmax < 0) = -9: The argument work[2] had an illegal value (iopt[7] = 1 and h0 = 0) = -10: The argument lwork had an illegal value (lwork too small) = -11: The argument iwork[0] had an illegal value (ml < 0 or ml > n) = -11: The argument iwork[1] had an illegal value (mu < 0 or mu > n) = -11: The argument iwork[2] had an illegal value (maxord not in range) = -12: The argument liwork had an illegal value (liwork too small) = -13: The argument info had an illegal value (info != 0, 1, 2 nor 11 to 13) = -17: The argument ldyypd had an illegal value (ldyypd < n and info(6) = 0, or ldyypd < 2*ml+mu+1 and info(6) = 1) = 1: Successful exit (t = tout). User can reenter with a new tout, which must be different from t. = 2: Interruption in intermediate-output mode (still not reached tout). User can reenter to resume for another step in the direction of tout. = 11: Maximum number of steps (500) exceeded. User can reenter to resume. Additional 500 steps will be allowed. = 12: Error tolerances are too stringent. User can reenter to continue with automatically relaxed tolerances. User can also reenter with manually changed tolerances. = 13: Pure relative error test (atol = 0) is impossible because computed solution is zero. User can reenter with positive atol value to continue computation. = 14: Repeated error test failure. = 15: Repeated corrector convergence test failure. = 16: The matrix of partial derivatives is singular. = 17: Multiple convergence test failure in this step. = 18: The corrector could not converge because ires = -1. = 19: ires = -2 was encountered. = 20: Failed to compute the initial yprime. = 21: Infinite loop has been detected. [out] tt irev = 1 to 7: The value of t where the residual values should be evaluated and given in yyp[] in the next call. [in] yyp[] Array yyp[lyyp] (lyyp >= n) irev = 1 to 7: The computed residual values at given t (= tt), y (= y[]) and y' (= yp[]), i.e. yyp[i] = fi(t, y, y') (i = 0 to n-1), should be given in the next call. [in] ldyypd Leading dimension of the two dimensional array yypd[][]. (ldyypd >= n (when iopt[5] = 0), ldyypd >= 2*ml + mu + 1 (when iopt[5] = 1)) [in,out] yypd[][] Array yypd[lyypd][ldyypd] (lyypd >= n) irev = 30 to 33: Two dimensional partial derivatives should be stored in general matrix form (when iopt[5] = 0) or band matrix form (when iopt[5] = 1). See below for further details of band matrix form. Since all elements of yypd[][] are initialized to 0 before returning with irev = 30 to 33, only the non-zero elements can be stored. [out] cj irev = 30 to 33: A scalar which is used to compute the matrix of partial derivatives. [in,out] ires irev = 1 to 7: An integer flag which is always equal to 0 on return from dassl. User should alter ires only if user encounters an illegal value of y[] or a stop condition. Set ires = -1 if a value of y[] is illegal, and dassl will try to solve the problem without getting ires = -1. If ires = -2, dassl will return control to the calling program with info = 11. [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 7: User should set the computed residual values f(t, y, y') at tt, y[] and yp[] in yyp[]. Do not alter any variables other than yyp[]. = 30 to 33: User should set the computed matrix of partial derivatives df/dy + cj*df/dy' at tt, y[] and yp[] in yypd[][]. Do not alter any variables other than yypd[]. Further Details The special band matrix form with 2*ml+mu+1 rows x n columns used in this routine 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 * * The first ml rows marked + are used as the work area. The elements marked * are not used by the routine. Reference SLATEC (DASSL)