Scilab Home page | Wiki | Bug tracker | Forge | Mailing list archives | ATOMS | File exchange
Change language to: English - Português - 日本語

See the recommended documentation of this function

# dae

Differential algebraic equations solver

### Calling Sequence

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

### 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 `ydot0` is the initial state derivative value or an estimation of it (see below).

t0

a real number, the initial time.

t

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

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

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

res

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

a Scilab function

In this case, Its calling sequence 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`, `=-1` if residue is locally not defined for `(t,x,xdot)`, `=-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 calling sequence 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 calling sequence 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 calling sequence 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 calling sequence 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 calling sequence 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 calling sequence 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 calling sequence of the function `surface` is now

`r=surface(t,x,p1,p2,...)`
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```

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

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

`integer nx, ng,ipar(*)`

• The C calling sequence must be

C2F(<r_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]` `times` is the value of the time at which the surface is crossed, `num` is the number of the crossed surface

hd

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

y

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 had been computed. Else `y` is the vector `[x(t);xdot(t)]`.

### Description

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

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

### Examples

```//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 with C code (c compiler needed) --------------------------------------------------
//-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 = yd - y;'
' res = yd - (100.0*(1.0 - y*y)*y - y);}'
' '
'void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)'
'{pd=*cj - 0.0;'
' pd=    - (-200.0*y*y - 1.0);'
' pd=    - 1.0;'
' pd=*cj - (100.0*(1.0 - y*y));}'
' '
'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
'{ groot = y;}']
mputl(code,TMPDIR+'/t22.c')
//-3- run
rtol=[1.d-6;1.d-6];atol=[1.d-6;1.d-4];
t0=0;y0=[2;0];y0d=[0;-2];t=[20:20:200];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)))```