dae

Differential algebraic equations solver

Syntax

```y = dae(initial, t0, t, res)
[y [,hd]] = dae(initial, t0, t [[,rtol], atol], res [,jac] [,hd])
[y, rd] = dae("root", initial, t0, t, res, ng, surface)
[y, rd [,hd]] = dae("root", initial, t0, t [[,rtol], atol], res [,jac], ng, surface [,hd])
[y, rd] = dae("root2", initial, t0, t, res, ng, surface)
[y, rd [,hd]] = dae("root2", initial, t0, t [[,rtol], atol], res [,jac], ng, surface [, psol, pjac] [, hd])```

Arguments

initial

a column vector. It may be equal to `x0` or `[x0;xdot0]`. Where `x0` is the state value at initial time `t0` and `xdot0` is the initial state derivative value or an estimation of it (see below).

t0

a real number, the initial time.

t

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

atol

a real scalar or a column vector of same size as `x0`, the absolute error tolerance of solution. If `atol` is a vector the tolerances are specified for each component of the state.

rtol

a real scalar or a column vector of same size as `x0`, the relative error tolerance of solution. If `rtol` is a vector the tolerances are specified for each component of the state.

res

an external function, computes the value of `g(t,y,ydot)`. It may be

a Scilab function

In this case, the syntax must be `[r,ires]=res(t,x,xdot)` and `res` must return the residue `r=g(t,x,xdot)` and error flag `ires`.

`ires = 0` if `res` succeeds to compute `r`.

`ires = -1` if residue is locally not defined for `g(t,x,xdot)`.

`ires =-2` if parameters are out of admissible range.

a list

This form of external is used to pass parameters to the function. It must be as follows:

```list(res,p1,p2,...)
```

where the syntax of the function `res` is now

```r=res(t,y,ydot,p1,p2,...)
```

`res` still returns the residual value as a function of `(t,x,xdot,x1,x2,...)`, and `p1, p2,...` are function parameters.

a character string

it must refer to the name of a C or Fortran routine. Assuming that <`r_name`> is the given name,

• The Fortran calling sequence must be

`<r_name>(t,x,xdot,res,ires,rpar,ipar)`

```double precision t,x(*),xdot(*),res(*),rpar(*)```

`integer ires,ipar(*)`

• The C syntax must be

```C2F(<r_name>)(double *t, double *x, double *xdot, double *res, integer *ires, double *rpar, integer *ipar)```

where

• `t` is the current time value

• `x` the state array

• `xdot` the array of state derivatives

• `res` the array of residuals

• `ires` the execution indicator

• `rpar` is the array of floating point parameter values, needed but cannot be set by the `dae` function

• `ipar` is the array of floating integer parameter values, needed but cannot be set by the `dae` function

jac

an external, computes the value of `dg/dx+cj*dg/dxdot` for a given value of parameter `cj`. It may be

a Scilab function

Its syntax must be `r=jac(t,x,xdot,cj)` and the `jac` function must return `r=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot` where `cj` is a real scalar.

a list

This form of external is used to pass parameters to the function. It must be as follows:

```list(jac,p1,p2,...)
```

where the syntax of the function `jac` is now

```r = jac(t,x,xdot,p1,p2,...)
```

`jac` still returns `dg/dx+cj*dg/dxdot` as a function of `(t,x,xdot,cj,p1,p2,...)`.

a character string

it must refer to the name of a C or Fortran routine. Assuming that <j_name> is the given name,

• The Fortran calling sequence must be

```<j_name>(t, x, xdot, r, cj, ires, rpar, ipar)```

double precision ```t, x(*), xdot(*), r(*), ci, rpar(*)```

integer `ires, ipar(*)`

• The C syntax must be

```C2F(<j_name>)(double *t, double *x, double *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)```

where `t, x, xdot, ires, rpar, ipar` have similar definition as above, `r` is the results array

surface

an external, computes the value of the column vector `surface(t,x)` with `ng` components. Each component defines a surface.

a Scilab function

Its syntax must be `r=surface(t,x)`, this function must return a vector with `ng` elements.

a list

This form of external is used to pass parameters to the function. It must be as follows:

```list(surface,p1,p2,...)
```

where the syntax of the function `surface` is now

```r = surface(t,x,p1,p2,...)
```
a character string

it must refer to the name of a C or Fortran routine. Assuming that <s_name> is the given name,

• The Fortran calling sequence must be

```<s_name>(nx, t, x, ng, r, rpar, ipar)```

```double precision t, x(*), r(*), rpar(*)```

`integer nx, ng,ipar(*)`

• The C syntax must be

```C2F(<s_name>)(double *t, double *x, double *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)```

where `t, x, rpar, ipar` have similar definition as above, `ng` is the number of surfaces, `nx` the dimension of the state and `r` is the results array.

rd

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

pjac

external (function, list or string). See daskr

psol

external (function, list or string). See daskr

hd

a real vector, as an output it stores the `dae` context. It can be used as an input argument to resume integration (hot restart).

y

a real matrix. If `%DAEOPTIONS(2)=1` , each column is the vector `[t;x(t);xdot(t)]` where `t` is time index for which the solution has been computed. Else `y` is the vector `[x(t);xdot(t)]`.

Description

The `dae` function is a gateway built above the dassl, dasrt and daskr functions designed for implicit differential equations integration.

Option `"root"` calls the dasrt routine, and `"root2"` calls daskr.

```g(t, x, xdot) = 0
x(t0) = x0 and xdot(t0) = xdot0
```

If `xdot0` is not given in the initial argument, the `dae` function tries to compute it solving `g(t,x0,xdot0)=0`.

If `xdot0` is given in the initial argument it may be either a compatible derivative satisfying `g(t,x0,xdot0)=0` or an approximate value. In the latter case %DAEOPTIONS(7) must be set to 1.

Detailed examples using Scilab and C coded externals are given in `modules/differential_equations/tests/unit_tests/dassldasrt.tst` and `modules/differential_equations/tests/unit_tests/daskr.tst`

Examples

Example #1: dassl (no root finding)

```// Example with Scilab code
//--------------------------------------------------
function [r, ires]=chemres(t, y, yd)
r(1) = -0.04*y(1) + 1d4*y(2)*y(3) - yd(1);
r(2) =  0.04*y(1) - 1d4*y(2)*y(3) - 3d7*y(2)*y(2) - yd(2);
r(3) =       y(1) +     y(2)      + y(3)-1;
ires =  0;
endfunction

function pd=chemjac(x, y, yd, cj)
pd = [-0.04-cj , 1d4*y(3)               , 1d4*y(2);
0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
1       , 1                      , 1       ]
endfunction

x0 = [1; 0; 0];
xd0 = [-0.04; 0.04; 0];
t = [1.d-5:0.02:.4, 0.41:.1:4, 40, 400, 4000, 40000, 4d5, 4d6, 4d7, 4d8, 4d9, 4d10];

y = dae([x0, xd0], 0, t, chemres); // Returns requested observation time points

%DAEOPTIONS = list([], 1, [], [], [], 0, 0); // Ask  dae mesh points to be returned
y = dae([x0, xd0], 0, 4d10, chemres); // Without jacobian
y = dae([x0, xd0], 0, 4d10, chemres, chemjac); // With jacobian```

Example #2: dasrt ("root")

```// Example with C code (C compiler needed)
//--------------------------------------------------
bOK = haveacompiler();
if bOK <> %t
[btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info');
return
end

//-1- Create the C codes in TMPDIR - Vanderpol equation, implicit form
code = ['#include <math.h>'
'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
'{res[0] = yd[0] - y[1];'
' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
' '
'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
'{pd[0] = *cj - 0.0;'
' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
' pd[2] =     - 1.0;'
' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
' '
'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
'{ groot[0] = y[0];}']
previous_dir = pwd();
cd TMPDIR;
mputl(code, 't22.c')

//-3- Run
rtol = [1.d-6; 1.d-6];
atol = [1.d-6; 1.d-4];
t0 = 0; t = [20:20:200];
y0 = [2; 0]; y0d = [0; -2];
ng = 1;

// Simple simulation
t = 0:0.003:300;
yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
clf(); plot(yy(1, :), yy(2, :))
// Find first point where yy(1) = 0
[yy, nn, hotd] = dae("root", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22');
plot(yy(1, 1), yy(2, 1), 'r+')
xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));

// Hot restart for next point
t01 = nn(1);
[pp, qq] = size(yy);
y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
[yy, nn, hotd] = dae("root", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', hotd);
plot(yy(1, 1), yy(2, 1), 'r+')
xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
cd(previous_dir);```

Example #3: daskr ("root2"), using default 'psol' and 'pjac' routines

```// Example with C code (C compiler needed)
//--------------------------------------------------
bOK = haveacompiler();
if bOK <> %t
[btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info');
return
end

//-1- Create the C codes in TMPDIR - Vanderpol equation, implicit form
code = ['#include <math.h>'
'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
'{res[0] = yd[0] - y[1];'
' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
' '
'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
'{pd[0] = *cj - 0.0;'
' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
' pd[2] =     - 1.0;'
' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
' '
'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
'{ groot[0] = y[0];}']
previous_dir = pwd();
cd TMPDIR;
mputl(code, 't22.c')

//-3- Run
rtol = [1.d-6; 1.d-6];
atol = [1.d-6; 1.d-4];
t0 = 0; t = [20:20:200];
y0 = [2; 0]; y0d = [0; -2];
ng = 1;

// Simple simulation
t = 0:0.003:300;
yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
clf(); plot(yy(1, :), yy(2, :))
// Find first point where yy(1) = 0
%DAEOPTIONS = list([] , 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
[yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
plot(yy(1, 1), yy(2, 1), 'r+')
xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));

// Hot restart for next point
t01 = nn(1);
[pp, qq] = size(yy);
y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
[yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
plot(yy(1, 1), yy(2, 1), 'r+')
xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
cd(previous_dir);```