daskr

Arguments

x0

is either y0 (ydot0 is estimated by daskr with zero as first estimate) or the matrix [y0 ydot0]. g(t, y0, ydot0) must be equal to zero. If you only know an estimate of ydot0, set info(7)=1.

y0

a real column vector of initial conditions.

ydot0

a real column vector of the time derivative of y at t0 (may be an estimate).

t0

a real number, it is the initial instant.

t

a real scalar or vector. Gives instants for which you want the solution. Note that you can get solution at each daskr's step point by setting info(2)=1.

atol, rtol

real scalars or column vectors of same size as y or both of size 1. atol and rtol give respectively absolute and relative error tolerances of solution. If vectors, the tolerances are specified for each component of y.

res

external (function, list or string). Computes the value of g(t, y, ydot). It may be :

• A Scilab function.

  Its syntax must be [r, ires] = res(t, y, ydot) and must return the residual r = g(t, y, ydot) and error flag ires. ires = 0 if res succeeds to compute r, ires = -1 if residual is locally not defined for (t, y, ydot), ires = -2 if parameters are out of admissible range.

• A list.

  This form allows to pass parameters other than t, y, ydot to the function. It must be as follows:

  list(res, x1, x2, ...)
  

  where the syntax of the function res is now

  r = res(t, y, ydot, x1, x2, ...)
  

  res still returns r = g(t, y, ydot) as a function of (t, y, ydot, x1, x2, ...).

  Warning: this form must not be used if there is no extra argument to pass to the function.

• A string.

  It must refer to the name of a C function or a Fortran subroutine linked with Scilab.

  In C, the syntax must be:

  void res(double *t, double y[], double yd[], double r[],
           int *ires, double rpar[], int ipar[])
  

  In Fortran, it must be:

  subroutine res(t, y, yd, r, ires, rpar, ipar)
  double precision t, y(*), yd(*),r(*),rpar(*)
  integer ires, ipar(*)
  

  The rpar and ipar arrays must be present but cannot be used.

jac

external (function, list or string). Computes the value of dg/dy + cj*dg/dydot for a given value of parameter cj.

• A Scilab function.

  Its syntax must be r = jac(t, y, ydot, cj) and must return r = dg(t, y, ydot)/dy + cj*dg(t, y, ydot)/dydot where cj is a real scalar.

• A list.

  It must be as follows

  list(jac, x1, x2, ...)
  

  where the syntax of the function jac is now

  r = jac(t, y, ydot, cj, x1, x2,...)
  

  jac still returns dg/dy + cj*dg/dydot as a function of (t, y, ydot, cj, x1, x2, ...).

• A character string.

  It must refer to the name of a C function or a Fortran subroutine linked with Scilab

  In C, the syntax must be:

  void jac(double *t, double y[], double yd[], double pd[],
           double *cj, double rpar[], int ipar[])
  

  In Fortran, it must be:

  subroutine jac(t, y, yd, pd, cj, rpar, ipar)
  double precision t, y(*), yd(*), pd(*), cj, rpar(*)
  integer ipar(*)
  

surf

external (function, list or string). Computes the value of the column vector surf(t, y) with ng components. Each component defines a surface. It may be defined by:

• A Scilab function.

  Its syntax must be surf(t, y)

• A list.

  It must be as follows

  list(surf, x1, x2, ...)
  

  where the syntax of the function surf is now

  r = surf(t, y, x1, x2, ...)
  

• A character string.

  It must refer to the name of a C function or a Fortran subroutine linked with Scilab.

  In C, the syntax must be:

  void surf(int *ny, double *t, double y[], int *ng, double gout[])
  

  In Fortran, it must be:

  subroutine surf(ny, t, y, ng, gout)
  double precision t, y(*), gout(*)
  integer ny, ng
  

info

list which contains 14 elements. Default value is list([], 0, [], [], [], 0, [], 0, [], 0, 0, [], [], 1).

info(1)

real scalar which gives the maximum time for which g is allowed to be evaluated or an empty matrix [] if no limits imposed for time.

info(2)

flag which indicates if daskr returns its intermediate computed values (= 1) or only the user specified time point values (= 0).

info(3)

2 components vector which give the definition [ml,mu] of band matrix computed by jac; r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)". If jac returns a full matrix set info(3)=[]. Treat as dummy if info(8)=1.

info(4)

real scalar which gives the maximum step size. Set info(4)=[] if no limitation.

info(5)

real scalar which gives the initial step size. Set info(5)=[] if not specified.

info(6)

set info(6)=1 if the solution is known to be non negative, else set info(6)=0.

info(7)

if ydot0 is set so that g(t0, y0, ydot0) = 0 then set info(7)=[]. Otherwise, set info(7)=[+-1, ..., +-1], with info(7)(i) = 1 if y(i) is a differential variable and info(7)(i) = -1 if y(i) is an algebraic variable (if its derivatives do not appear explicitly in the function g(t, y, ydot)).

info(8)

direct / Krylov method. Set info(8)=1 and provide a routine in psol if you want the solver to use Krylov iterations, else (daskr's direct method) set info(8)=0.

info(9)

Krylov parameters. Treat as dummy argument if you have set info(8)=0. Otherwise, set info(9)=[] or info(9)=[maxl kmp nrmax epli], where:

- maxl = maximum number of iterations in the GMRes algorithm (default min(5, neq)),

- kmp = number of vectors on which orthogonalization is done in the GMRes algorithm (default maxl),

- nrmax = maximum number of restarts of the GMRes algorithm per nonlinear iteration (default 5),

- epli = convergence test constant in GMRes algorithm (default 0.05).

info(10)

initial conditions. Treat as dummy argument if info(7)=[]. Set info(10)=1 if the solver should stop right after computation of the initial conditions, else set info(10)=0.

info(11)

preconditioner computation and LU-factorization routine for psol. Treat as dummy argument if info(8)=0. Set info(11)=1 and provide a pjac routine if the external psol should use a specific routine, else set info(11)=0.

info(12)

if you wish to control errors locally on all the variables then set info(12)=[]. Otherwise, set info(12)=[+-1, ..., +-1], with info(12)(i) = 1 if y(i) is a differential variable and info(12)(i) = -1 if y(i) is an algebraic variable (if its derivatives do not appear explicitly in the function g(t, y, ydot)).

info(13)

heuristic parameters. Treat as dummy argument if info(7)=[]. Otherwise, set info(13)=[] or info(13)=[mxnit mxnj mxnh lsoff stptol epinit], where:

- mxnit = maximum number of Newton iterations per Jacobian or preconditioner evaluation (default 5 if info(8)=0, 15 otherwise),

- mxnj = maximum number of Jacobian or preconditioner evaluations (default 6 if info(8)=0, 2 otherwise),

- mxnh = maximum number of values of the artificial stepsize parameter h to be tried if info(7) ≠ [] (default 5),

- lsoff = flag to turn off the linesearch algorithm (lsoff = 0 means linesearch is on, lsoff = 1 means it is turned off) (default 0),

- stptol = minimum scaled step in linesearch algorithm (default (unit roundoff)^(2/3)),

- epinit = swing factor in the Newton iteration convergence test (default 0.01).

info(14)

verbosity. Set info(14)=1 for minimal extra printing, info(14)=2 for full printing, else set info(14)=0.

psol

external (function, list or string). Solves a linear system P*x = b, with P being the factored preconditioner that routine pjac computed beforehand and stored in wp and iwp.

• A Scilab function.

  Its syntax must be [r, ier] = psol(wp, iwp, b) and must return the solution of the system in r and an error flag ier.

• A list.

  It must be as follows:

  list(psol, x1, x2, ...)
  

  where the syntax of psol is now

  psol(wp, iwp, b, x1, x2, ...)
  

  psol still returns the solution in r.

• A character string.

  It must refer to the name of a C function or a Fortran subroutine linked with Scilab

  In C, the syntax must be:

  void psol (int*neq, double*t, double*y, double*ydot, double*savr,
    double*wk, double*cj, double*wght, double*wp, int*iwp, double*b, double*eplin, int*ier, double*rpar, int*ipar)
  

  where the arrays wp and iwp contain matrix elements of LU-factored preconditioner P, wp being the values and iwp the pivots used in the factorization.

  In Fortran, it must be:

  subroutine psol (neq, t, y, ydot, savr, wk, cj, wght,
    wp, iwp, b, eplin, ier, rpar, ipar)
  double precision t,y(*), ydot(*), savr(*), wk(*), cj, wght(*), wp(*),
    b(*), eplin, rpar(*)
  integer neq, iwp(*), ier, ipar(*)
  

pjac

external (function, list or string). Computes the value of dg/dy + cj*dg/dydot for a given value of parameter cj and LU-factorizes it in two arrays, real and integer.

• A Scilab function.

  Its syntax must be [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr) and in return, the arrays wp and iwp must contain all factored preconditioner information.

• A list.

  It must be as follows

  list(pjac, x1, x2, ...)
  

  where the syntax of pjac is

  pjac(neq, t, y, ydot, h, cj, rewt, savr, x1, x2,...)
  

  pjac still returns factorized dg/dy + cj*dg/dydot as a function of (neq, t, y, ydot, h, cj, rewt, savr, x1, x2, ...).

• A character string.

  It must refer to the name of a C function or a Fortran subroutine linked with Scilab

  In C, the syntax must be:

  void pjac (double*res, int*ires, int*neq, double*t, double*y, double*ydot, double*rewt, double*savr,
  double*wk, double*h, double*cj, double*wp, int*iwp, int*ier, double*rpar, int*ipar)
  

  In Fortran, it must be:

  subroutine pjac (res, ires, neq, t, y, ydot, rewt, savr,
    wk, h, cj, wp, iwp, ier, rpar, ipar)
  double precision res(*), t, y(*), ydot(*), rewt(*), savr(*),
    wk(*), h, cj, wp(*), rpar(*)
  integer ires, neq, iwp(*), ier, ipar(*)
  

hd

real vector which allows to store the daskr context and to resume integration.

r

real matrix. Each column is the vector [t; x(t); xdot(t)] where t is the time index for which the solution has been computed, x(t) is the value of the computed solution, xdot(t) is the derivative of the computed solution.

nn

a vector with two entries [times num], times is the value of the time at which the surface is crossed, num is the number of the crossed surface.

Description

Solution of the implicit differential equation:

g(t, y, ydot) = 0
y(t0) = y0  and  ydot(t0) = ydot0

Returns the surface crossing instants and the number of the surface reached in nn.

Detailed examples can be found in SCI/modules/differential_equations/tests/unit_tests/daskr.tst

Examples

// dy/dt = ((2*log(y)+8)/t-5)*y,  y(1) = 1,  1 <= t <= 6
// g1 = ((2*log(y)+8)/t-5)*y
// g2 = log(y) - 2.2491
y0 = 1; t = 2:6; t0 = 1; y0d = 3;
atol = 1.d-6; rtol = 0; ng = 2;

deff('[delta, ires] = res1(t, y, ydot)', 'ires = 0; delta = ydot-((2*log(y)+8)/t-5)*y')
deff('[rts] = gr1(t, y)', 'rts = [((2*log(y)+8)/t-5)*y; log(y)-2.2491]')

[yy, nn] = daskr([y0, y0d], t0, t, atol, rtol, res1, ng, gr1);
nn

// Should return nn = [2.4698972 2]

See also

• ode — ordinary differential equation solver
• dasrt — DAE solver with zero crossing
• dassl — differential algebraic equation
• impl — differential algebraic equation
• call — Fortran or C user routines call
• link — dynamic linker
• external — Scilab Object, external function or routine

History