XLPack 7.0
XLPack Numerical Library (C API) Reference Manual
Loading...
Searching...
No Matches
debdf
debdf_r
radaua
radaua_r
rodasa
rodasa_r
seulexa
seulexa_r

◆ radaua()

void radaua ( int  n,
void(*)(int, double, double *, double *)  f,
double *  t,
double  y[],
double  tout,
double  tend,
double *  rtol,
double *  atol,
int  itol,
void(*)(int, double, double *, int, double *)  fjac,
int  ijac,
int  mljac,
int  mujac,
void(*)(int, int, double *)  fmas,
int  imas,
int  mlmas,
int  mumas,
int  mode,
double  work[],
int  lwork,
int  iwork[],
int  liwork,
int *  info 
)

Initial value problem of ordinary differential equations (variable (5, 9, 13-th) order implicit Runge-Kutta method (Radau IIA))

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 variable (5, 9, 13-th) order implicit Runge-Kutta method (Radau IIA)) code, RADAU (Reference (1)).
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[])
{
Compute dy/dt and store to yp[].
}
where n is the number of equations, and dy/dt (= f(t, y[])) is the computed derivative at given t and y[].
[in,out]tIndependent 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]toutSet 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]tendThe 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]rtolScalar 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]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) 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]itolSpecifies 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]fjacThe user supplied subroutine, which computes the Jacobian, defined as follows.
void fjac(int n, double t, double y[], int ldypd, double ypd[][ldypd])
{
Compute Jacobian at t and y[] and store in ypd[][].
}
Jacobian should be stored in ypd[][] in full matrix form (general n x n two dimensional array) if mljac = n, or in band matrix form (see details below) if mljac < n.
[in]ijacSpecify how to compute the Jacobian.
= 0: Jacobian is computed by finite difference approximation. fjac is not referenced.
= 1: Jacobian is computed by flac.
[in]mljacThe 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]mujacThe upper bandwidth of Jacobian. (0 <= mujac <= n)
If mljac = n, mujac is ignored.
[in]fmasThe user supplied subroutine, which provides the mass matrix M, defined as follows.
void fmas(int n, int ldam, double am[][ldam])
{
Store mass matrix in am[][].
}
Mass matrix should be stored in am[][] in full matrix form if mlmas = n, or in band matrix form if mlmas < n.
[in]imasSpecify whether mass matrix M is identity matrix.
= 0: M is the identity matrix. fmas is not referenced.
= 1: M is supplied by fmas.
[in]mlmasThe 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]mumasThe upper bandwidth of mass matrix M. (0 <= mumas <= n)
If mlmas = n, mumas is ignored.
[in]modeMode 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 (e.g. 0.001). Negative value forces to compute the Jacobian after every accepted step.
work[4] and work[5]: Parameters for step size selection. (work[4] <= 1, work[5] >= 1. (default work[4] = 0.2, work[5] = 8)
The new step size is chosen subject to the restriction work[4] <= hnew/hold <= work[5].
work[8]: The safety factor in step size prediction (0.001 < work[8] < 1). (default = 0.9)
work[11], work[12]: If work[11] < hnew/hold < work[12], the step size is not changed. (work[11] <= 1, work[12] >= 1) (default: work[11] = 1, work[12] = 1.2)
For large systems, work[11] = 0.99 and work[12] = 2 might be good. This saves, together with a large work[2], LU decompositions and computing time. For small systems default values might be good.
work[13]: Order is increased if the contractivity factor is larger than this value. (default = 0.002)
work[14]: Order is decreased if the contractivity factor is larger than this value. (default = 0.8)
work[15], work[16]: Order is decreased only if the step size ratio safisfies work[15] <= hnew/h <= work[16]. (default: work[15] = 1.2, work[16] = 0.8)
[out]
work[1]: Last step size.
[in]lworkSize of array work[]. (abs(lwork) >= n*(3 + 3*nsmax + ldjac) + nm1*(nsmax*lde + ldmas) + 3*nsmax + 60 + n*(max(ldjac, ldmas) + 2))
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])
nsmax = iwork[10]/10
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 = 10000)
iwork[2]: Maximum number of Newton iterations for the solution of the implicit system in each step. (default = 7)
iwork[3]: Specifies the starting value of Newton iteration. (default = 0)
= 0: Use extrapolated solution as starting value.
= 1: Use 0 as starting value.
The latter is recommended if Newton method has difficulties with convergence.
iwork[4], iwork[5], iwork[6]: Dimensions of the index 1, 2 and 3 variables. (iwork[4] > 0, iwork[4] + iwork[5] + iwork[6] = n) (default: iwork[4] = n, iwork[5] = 0, iwork[6] = 0)
These parameters are important for differential algebraic equations (DAEs) of index > 1. The function subroutine should be written such that the index 1, 2 and 3 variables appear in this order. In estimating the error, the index 2 variables are multiplied by h, the index 3 variables by h^2. For ordinary defferential equations (ODEs), iwork[4] equals to n.
iwork[7]: Specifies the step size strategy. (default = 1)
= 1: Model predictive controller (Gustafsson).
= 2: Classical step size control.
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[10]: Minimal number of stages nsmin and maximal number of stages nsmax (= 3, 5 or 7). 10*nsmax + nsmin should be set. (default = 73 (nsmin = 3, nsmax = 7))
ns correspond to implicit Runge-Kutta method of order 2*ns - 1.
iwork[11]: Value of ns for the first step. (3, 5 or 7) (default = nsmin)
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]liworkSize of array iwork[]. (abs(liwork) >= (2 + (nsmax - 1)/2)*nm1 + 80, nsmax = iwork[10]/10, nm1 = n - m1 (m1 = iwork[8]))
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.
= 12: (Error) Step size becomes too small.
= 13: (Error) Matrix is repeatedly singular.
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)
Generated by doxygen 1.9.7