differential algebraic equation
[r [,hd]]=dassl(x0,t0,t [,atol,[rtol]],res [,jac] [,info] [,hd])
: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 a real column vector of initial conditions. : :ydot0 a real column vector of the time derivative of y at t0
(may be an estimate).:
: :t0 a real number is the initial instant. : :t a 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.
: :res an 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 an 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 C or 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) a 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) a flag which indicates if dassl returns its intermediate
- computed values ( flag=1) or only the user specified time point values ( flag=0).
- : :info(3) a 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) a real scalar which gives the maximum step size. Set
- info(4)=[] if no limitation.
- : :info(5) a real scalar which gives the initial step size. Set
- info(5)=[] if not specified.
- : :info(6) set info(6)=1 if the solution is known to be non
- negative, else set info(6)=0.
- : :info(7) set info(7)=1 if ydot0 is just an estimation,
- info(7)=0 if g(t0,y0,ydot0)=0.
:
:
The dassl function integrate the differential algebraic 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;
y1=dassl([y0,yd0],0,4d10,chemres,info);
y2=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
y3=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],`fullfile`_(TMPDIR,'mycode.c')) //create the C file
// - compile it
`ilib_for_link`_(['chemres','chemjac'],`fullfile`_(TMPDIR,'mycode.c'),[],'c','',`fullfile`_(TMPDIR,'loader.sce'));//compile
// - link it with Scilab
`exec`_(`fullfile`_(TMPDIR,'loader.sce')) //incremental linking
// - call dassl
y4=dassl([y0,yd0],0,t,'chemres','chemjac');