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

See the recommended documentation of this function

# dassl

differential algebraic equation

### Calling Sequence

`[r [,hd]]=dassl(x0,t0,t [,atol,[rtol]],res [,jac] [,info] [,hd])`

### Arguments

x0

is either `y0` (`ydot0` is estimated by `dassl` 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

real column vector of initial conditions.

ydot0

real column vector of the time derivative of `y` at `t0` (may be an estimate).

t0

real number is the initial instant.

t

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

atol,rtol

real scalars or column vectors of same size as `y`. `atol,rtol` give respectively absolute and relative error tolerances of solution. If vectors the tolerances are specified for each component of `y`.

res

external (function or list or string). Computes the value of `g(t,y,ydot)`. It may be :

• A Scilab function.

Its calling sequence must be `[r,ires]=res(t,y,ydot)` and `res` must return the residue `r=g(t,y,ydot)` and error flag `ires`. `ires = 0` if `res` succeeds to compute `r`, `=-1` if residue is locally not defined for `(t,y,ydot)`, `=-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 calling sequence 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,...)`.

• A string.

it must refer to the name of a C or fortran subroutine linked with Scilab.

In C The calling sequence 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 or list or string). Computes the value of `dg/dy+cj*dg/dydot` for a given value of parameter `cj`

• A Scilab function.

Its calling sequence must be `r=jac(t,y,ydot,cj)` and the `jac` function 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 calling sequence 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 fortran subroutine linked with scilab

In C The calling sequence 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(*)```
info

optional list which contains `7` elements. Default value is list([],0,[],[],[],0,0);

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 `dassl` returns its intermediate computed values (`flag=1`) or only the user specified time point values (`flag=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)=[]`.

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(4)=[]` if not specified.

info(6)

set `info(6)=1` if the solution is known to be non negative, else set .

info(7)

set `info(7)=1` if `ydot0` is just an estimation, `info(7)=0` if `g(t0,y0,ydot0)=0`.

hd

real vector which allows to store the `dassl` context and to resume integration

r

real matrix . Each column is the vector [t;x(t);xdot(t)] where t is time index for which the solution had been computed

### Description

The dassl function integrate the algebro-differencial equation and returns the evolution of`y`a given time points

```g(t,y,ydot)=0
y(t0)=y0  and   ydot(t0)=ydot0```

### Examples

```function [r, ires]=chemres(t, y, yd)
r=[-0.04*y(1)+1d4*y(2)*y(3)-yd(1)
0.04*y(1)-1d4*y(2)*y(3)-3d7*y(2)*y(2)-yd(2)
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

y0=[1;0;0];
yd0=[-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=dassl([y0,yd0],0,t,chemres);

info=list([],0,[],[],[],0,0);
info(2)=1;
y=dassl([y0,yd0],0,4d10,chemres,info);
y=dassl([y0,yd0],0,4d10,chemres,chemjac,info);

//Using extra argument for parameters
//-----------------------------------
function [r, ires]=chemres(t, y, yd, a, b, c)
r=[-a*y(1)+b*y(2)*y(3)-yd(1)
a*y(1)-b*y(2)*y(3)-c*y(2)*y(2)-yd(2)
y(1)+y(2)+y(3)-1];
ires=0
endfunction

function pd=chemjac(x, y, yd, cj, a, b, c)
pd=[-a-cj , b*y(3)             , b*y(2);
a    ,-b*y(3)-2*c*y(2)-cj ,-b*y(2);
1    , 1                  , 1       ]
endfunction

y=dassl([y0,yd0],0,t,list(chemres,0.04,1d4,3d7),list(chemjac,0.04,1d4,3d7));

//using C code
//------------
// - create the C code
rescode=['void chemres(double *t, double y[], double yd[], double r[], int *ires, double rpar[], int ipar[])'
' {'
'   r[0] = -0.04*y[0]+1.0e4*y[1]*y[2]                -yd[0];'
'   r[1] =  0.04*y[0]-1.0e4*y[1]*y[2]-3.0e7*y[1]*y[1]-yd[1];'
'   r[2] =       y[0]+y[1]+y[2]-1;'
'   *ires = 0;'
' }'];

jaccode=['void chemjac(double *t, double y[], double yd[], double pd[], double *cj, double rpar[], int ipar[])'
' {'
'   /* first column*/'
'   pd[0] = -0.04-*cj;'
'   pd[1] =  0.04;'
'   pd[2] =  1.0;'
'    /* second column*/'
'   pd[3] =  1.0e4*y[2];'
'   pd[4] = -1.0e4*y[2]-2*3.0e7*y[1]-*cj;'
'   pd[5] =  1.0;'
'    /* third column*/'
'   pd[6] =  1.0e4*y[1];'
'   pd[7] = -1.0e4*y[1];'
'   pd[8] =  1.0;'
' }'];
mputl([rescode;jaccode],TMPDIR+'/mycode.c') //create the C file

// - compile it

// - link it with Scilab

// - call dassl
y=dassl([y0,yd0],0,t,'chemres','chemjac');```

• ode — ordinary differential equation solver
• dasrt — DAE solver with zero crossing
• impl — differential algebraic equation
• fort — Fortran or C user routines call