XLPack 6.1
C/C++ API Reference Manual
Loading...
Searching...
No Matches
_dassl
_dassl_r
_debdf
_debdf_r
_radau
_radau5
_radau5_r
_radau_r
_radaup
_radaup_r
_rodas
_rodas_r
_seulex
_seulex_r

◆ _rodas()

void _rodas ( int  n,
void(*)(int, double, double *, double *)  f,
int  ifcn,
double *  t,
double  y[],
double  tout,
double  rtol[],
double  atol[],
int  itol,
void(*)(int, double, double *, int, double *)  jac,
int  ijac,
int  mljac,
int  mujac,
void(*)(int, double, double *, double *)  dfx,
int  idfx,
void(*)(int, int, double *)  mas,
int  imas,
int  mlmas,
int  mumas,
void(*)(int, double, double, double *, int, double *, int *)  solout,
int  iout,
double  work[],
int  lwork,
int  iwork[],
int  liwork,
double  cont[],
int  lcont,
int *  info 
)

Initial value problem of ordinary differential equations (4(3)-th order Rosenbrock method)

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).

rodas is the code based on the 4(3)-th order Rosenbrock method. It is provided with the step control algorithm and the continuous output feature.

See for details in the reference below.
Parameters
[in]nNumber of differential equations. (n >= 1)
[in]fThe 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 many cases, M is an identity matrix. However, it can be defined in subroutine mas if necessary.
[in]ifcnGives 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]tThis 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]toutSet 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]rtolScalar 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]atolScalar 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]itolSpecifies 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]jacThe user supplied subroutine, which computes the Jacobian dfi(t, y)/dyj, defined as follows.
void jac(int n, double t, double y[], int ldypd, double ypd[][ldypd])
{
Store Jacobian at t and y[] in ypd[][].
}
Two dimensional Jacobian data should be stored in ypd[][] with the row dimension ldypd in general matrix form if mljac = n, or band matrix form if mljac < n.
See below for further details of band matrix form.
[in]ijacSwitch for the computation of the Jacobian.
= 0: Jacobian is computed internally by finite differences, jac is never called.
= 1: Jacobian is supplied by jac.
(For other values, ijac = 0 is assumed.)
[in]mljacThe 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]mujacThe upper bandwidth of Jacobian. (0 <= mujac <= n)
If mljac = n, mujac is ignored.
[in]dfxThe user supplied subroutine, which computes the partial derivatives of f(t, y) with respect to t, defined as follows.
void dfx(int n, double t, double y[], double fx[])
{
fx[i] = dfi/dt (for i = 0 to n-1)
}
This routine will be called only if idfx = 1 and ifcn = 1.
[in]idfxSwitch for the computation of the partial derivatives of f(t, y).
= 0: The partial derivatives are computed internally by finite defferences.
= 1: The partial derivatives are supplied by dfx.
(For other values, idfx = 0 is assumed.)
[in]masThe user supplied subroutine, which provides the mass matrix M, defined as follows.
void mas(int n, int ldam, double am[][ldam])
{
Store mass matrix in am[][].
}
Two dimensional mass matrix data should be stored in am[][] with the row dimension ldam in general matrix form if mlmas = n, or band matrix form if mlmas < n.
See below for further details of band matrix form.
[in]imasSwitch for the mass matrix M.
= 0: M is supposed to be the identity matrix. mas is never called.
= 1: M is supplied by mas.
(For other values, imas = 0 is assumed.)
[in]mlmasThe 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]mumasThe upper bandwidth of mass matrix M. (0 <= mumas <= n)
If mlmas = n, mumas is ignored.
[in]soloutThe user supplied subroutine to print out the intermediate solutions, which is called after every successful step, defined as follows.
void solout(int nr, double told, double t, double y[], int n, double cont[], int *irtrn)
{
Output the y[] values at nr-th step t. told is the previous value of t. n is the order of equations.
The input value of irtrn will be 0, 1 or 2 in the first, intermediate or last call of solout, respectively. irtrn also serves to interrupt the integration. If irtrn is set to the negative value in solout, the integration will be interrupted and exit with info = 2.
Dense output is supported by the control information cont[] if iout = 1.
The solution y[i] (0 <= i <= n-1) at the arbitrary point t2 in the interval [told, t] can be computed by the function call
y[i] = contro(i, t2, cont);
}
[in]ioutSwitch for calling the subroutine solout.
= 0: solout is never called
= 1: solout is used for output
(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]: Maximal step size. (default = tout - t)
work[3] and work[4]: Parameters for step size selection. (work[3] <= 1, work[4] >= 1)
The new step size is chosen subject to the restriction work[3] <= hnew/hold <= work[4]. (default work[3] = 0.2, work[4] = 6)
work[7]: The safety factor in step size prediction. (0.001 < work[7] < 1) (default = 0.9)
[out]
work[0]: Predicted step size of the last accepted step.
[in]lworkSize of array work[]. (lwork >= n*(ljac + ldypd + 13) + nm1*(lmas + le) + 20)
where
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])
ldypd = max(ljac, lmas)
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]: Switch for the choice of the coefficients. (default = 1)
= 1: Method in the reference page 452
= 2: Same method with different parameters
= 3: Method with coefficients of Gerd Steinebach
(For other values, default value will be used)
iwork[1]: Maximum number of allowed steps. (default = 100000) (if iwork[1] < 0, default value will be used)
iwork[7]: Switch for step size strategy. (default = 1)
= 1: Model predictive controller (Gustafsson)
= 2: Classical step size control
(For other values, default value will be used)
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]liworkSize of array iwork[]. (liwork >= nm1 + 20, where nm1 = n - m1 (m1 = iwork[8]))
If liwork < 0, abs(liwork) will be used and iwork[0] to iwork[19] will be initialized to zeros.
[out]cont[]Array cont[lcont]
Control information area for dense output. Used only when iout = 1.
[in]lcontSize of array cont[]. (lcont >= 4*n)
[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)
= -7: The argument rtol had an illegal value (rtol < 0 or rtol[i] < 0)
= -7: The argument rtol and atol had an illegal value (rtol = 0 and atol = 0, or rtol[i] = 0 and atol[i] = 0)
= -8: The argument atol had an illegal value (atol < 0 or atol[i] < 0)
= -18: The argument mljac or mlmas had an illegal value (mlmas > mljac)
= -19: The argument mujac or mumas had an illegal value (mumas > mujac)
= -22: The argument work[3] or work[4] had an illegal value (work[3] > 1 or work[4] < 1)
= -22: The argument work[7] had an illegal value (work[7] <= 0.001 or work[7] >= 1)
= -23: The argument lwork had an illegal value (lwork too small)
= -24: The argument iwork[0] had an illegal value (iwork[0] < 0 or iwork[0] > 3)
= -24: The argument iwork[8] or iwork[9] had an illegal value (iwork[8] < 0 or iwork[9] < 0 or iwork[8] + iwork[9] > n)
= -25: The argument liwork had an illegal value (liwork too small)
= -27: The argument lcont had an illegal value (lcont too small)
= -28: the argument info had an illegal value (info != 0 and info != 1)
= 1: Successful exit
= 2: interrupted by solout (normal return)
= 11: maximum number of steps exceeded
= 12: step size becomes too small
= 13: matrix is repeatedly singular
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[3] and iwork[4]).

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)
Generated by doxygen 1.9.6