# daskr

DAE solver with zero crossing

### Syntax

[r, nn [, hd]] = daskr(x0, t0, t [, atol [, rtol]], res [, jac], ng, surf [, info [, psol] [, pjac]] [, hd])

### 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 nowr = 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 nowr = 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 nowr = 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 nowpsol(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:

where the arraysvoid 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)

`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`

ispjac(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

### History

Version | Description |

5.5.0 | daskr solver added |

Report an issue | ||

<< daeoptions | Differential calculus, Integration | dasrt >> |