debdf debdf_r radaua radaua_r rodasa rodasa_r seulexa seulexa_r

debdf() void debdf ( int  n, void(*)(int, double, double *, double *)  f, double *  t, double  y[], double  tout, double *  rtol, double *  atol, int  itol, int  mode, int  itstop, void(*)(int, double, double *, int, double *)  djac, int  idjac, double  work[], int  lwork, int  iwork[], int  liwork, int *  info  ) Initial value problem of ordinary differential equations (1~5-th order backward differentiation formula (BDF)) 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. This is the code in the package of differential equation solvers DEPAC, consisting of the codes derkf, deabm, and debdf. debdf is a variable order (one through five) backward differentiation formula code. It is the most complicated of the three choices. debdf is primarily designed to solve stiff differential equations at crude to moderate tolerances. If the problem is very stiff at all, derkf and deabm will be quite inefficient compared to debdf. However, debdf will be inefficient compared to derkf and deabm on non-stiff problems because it uses much more storage, has a much larger overhead, and the low order formulas will not give high accuracies efficiently. debdf is a driver for a modification of the code LSODE written by A. C. Hindmarsh. 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. 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 mode. [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] tout Set tout to the point at which a solution is desired. You can take tout = t, in which case the code will evaluate the derivative of the solution at t and return. 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] 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 itol. 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 itol. 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] itol 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. (For other values, itol = 0 is assumed.) [in] mode Mode of operation. = 0: Return only at tout (interval mode) = 1: Return at every step (intermediate output mode) (For other values, mode = 0 is assumed.) [in] itstop 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 (For other values, itstop = 0 is assumed.) [in] djac The user supplied subroutine, which computes the Jacobian analytically, defined as follows. void djac(int n, double t, double y[], int ldypd, double ypd[][ldypd]) { dfi/dyj to be stored in ypd[][] } Two dimensional Jacobian should be stored in ypd[][] in general matrix form (when idjac = 0) or band matrix form (when idjac = 1). See below for further details of band matrix form. Since all elements of ypd[][] are initialized to 0 before calling djac, only the non-zero elements can be stored. If djac = NULL is specified, Jacobian is approximated by numerical differencing internally. The solution may be more reliable if the Jacobian is compute analytically by djac. Sometimes numerical differencing is cheaper than evaluating derivatives in djac analytically and sometimes it is not - this depends on the problem. [in] idjac Specify whether Jacobian is computed as full (general) matrix or band matrix. In general, computation with band matrix costs less time and storage than with full (general) 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) (For other values, idjac = 0 is assumed.) [out] work[] Array work[lwork] Work array. [in] lwork The length of work[]. (lwork >= 2*n^2 + 12*n + 24 (idjac = 0 (full Jacobian)), (3*ml + 2*mu + 14)*n + 24 (idjac = 1 (banded Jacobian))) [out] iwork[] Array iwork[liwork] Work array. [in] liwork The length of iwork[]. (liwork >= n + 38) [in,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 15: When returned from the routine with info = 1, 2 or 11 to 15, 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 15 as a next action if necessary. = -1: The argument n had an illegal value (n < 1) = -3: The argument t had an illegal value (t = tout or t != previous tout on continuation call) = -5: The argument tout had an illegal value (direction of integration is changed, or, tout conflicts with tstop) = -6: The argument rtol had an illegal value (rtol < 0 or rtol[i] < 0) = -7: The argument atol had an illegal value (atol < 0 or atol[i] < 0) = -14: The argument lwork had an illegal value (lwork too small) = -15: The argument iwork had an illegal value (ml < 0 or ml >= n, or mu < 0 or mu >= n) = -16: The argument liwork had an illegal value (liwork to small) = -17: The argument info had an illegal value (info != 0, 1, 2 nor 11 to 15) = 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 (10000) exceeded. User can reenter to resume. Additional 10000 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 convergence test failures on the last attempted step. User may reenter without changing info to continue, but an inaccurate Jacobian may be the problem. = 15: Repeated error test failures on the last attempted step. User may reenter without changing info to continue, but a singularity in the solution may be present. = 16: Infinite loop has been detected. 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. Reference SLATEC (DEPAC)