Scilab Website | Contribute with GitLab | Mailing list archives | ATOMS toolboxes
Scilab Online Help
2023.0.0 - Português

# 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 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 — Solucionador de equações diferenciais ordinárias
• dasrt — Solucionador de equações diferenciais algébricas com cruzamento de zeros
• dassl — Equação diferencial algébrica
• impl — equações diferenciais algébricas
• call — chamada a rotinas de usuário FORTRAN ou C
• link — dynamic linker
• external — objeto Scilab, função ou rotina externa

### History

 Versão Descrição 5.5.0 daskr solver added
 Report an issue << daeoptions Equações Diferenciais dasrt >>

 Copyright (c) 2022-2023 (Dassault Systèmes)Copyright (c) 2017-2022 (ESI Group)Copyright (c) 2011-2017 (Scilab Enterprises)Copyright (c) 1989-2012 (INRIA)Copyright (c) 1989-2007 (ENPC)with contributors Last updated:Tue Mar 07 09:28:40 CET 2023