dassl — differential algebraic equation
[r [,hd]]=dassl(x0,t0,t [,atol,[rtol]],res [,jac] [,info] [,hd])
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
real column vector of initial conditions.
real column vector of the time derivative of
y at t0 (may be an
estimate).
real number is the initial instant.
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.
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.
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.
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(*)
optional list which contains 7 elements.
Default value is list([],0,[],[],[],0,0);
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.
flag which indicates if dassl returns
its intermediate computed values (flag=1)
or only the user specified time point values
(flag=0).
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)=[].
real scalar which gives the maximum step size. Set
info(4)=[] if no limitation.
real scalar which gives the initial step size. Set
info(4)=[] if not specified.
set info(6)=1 if the solution is
known to be non negative, else set .
set info(7)=1 if
ydot0 is just an estimation,
info(7)=0 if
g(t0,y0,ydot0)=0.
real vector which allows to store the dassl
context and to resume integration
real matrix . Each column is the vector [t;x(t);xdot(t)] where t is time index for which the solution had been computed
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
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
ilib_for_link(['chemres','chemjac'],'mycode.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/loader.sce');//compile
// - link it with Scilab
exec(TMPDIR+'/loader.sce') //incremental linking
// - call dassl
y=dassl([y0,yd0],0,t,'chemres','chemjac');