Change language to:
Français - 日本語 - Português - Русский

Please note that the recommended version of Scilab is 2024.1.0. This page might be outdated.
However, this page did not exist in the previous stable version.

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

• 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