dae — Differential algebraic equations solver
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])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).
a real number, the initial time.
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.
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.
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.
an external. Computes the value of
g(t,y,ydot). It may be
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.
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.
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
an external. Computes the value of
dg/dx+cj*dg/dxdot for a given value of parameter
cj. It may be
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
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,...).
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
an external. Computes the value of the column vector
surface(t,x) with ng
components. Each component defines a surface.
Its calling sequence must be
r=surface(t,x), this function must return a
vector with ng elements.
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,...)
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.
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
a real vector, an an output it stores the
dae context. It can be used as an input argument
to resume integration (hot restart).
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)].
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
//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[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];}']
mputl(code,TMPDIR+'/t22.c')
//-2- compile and load them
ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/t22loader.sce');
exec(TMPDIR+'/t22loader.sce')
//-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)))