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

Please note that the recommended version of Scilab is 6.0.2. This page might be outdated.
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 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 .

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

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 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,...)
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 calling sequence 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

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

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;}']
cd TMPDIR;
mputl(code, 't22.c')
//-2- compile and load them
//-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)))

• ode — ordinary differential equation solver
• daeoptions — set options for dae solver
• dassl — differential algebraic equation
• impl — differential algebraic equation
• fort — Fortran or C user routines call