Please note that the recommended version of Scilab is 2024.1.0. This page might be outdated.
However, this page did not exist in the previous stable version.
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 bydassl
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 ofydot0
setinfo(7)=1
.- y0
a real column vector of initial conditions.
- ydot0
a real column vector of the time derivative of
y
att0
(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
.- atol, rtol
real scalars or column vectors of same size as
y
or both of size 1.atol, rtol
give respectively absolute and relative error tolerances of solution. If vectors, the tolerances are specified for each component ofy
.- 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)
andres
must return the residuer=g(t,y,ydot)
and error flagires
.ires = 0
ifres
succeeds to computer
,=-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 nowr=res(t,y,ydot,x1,x2,...)
res
still returnsr=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:
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
andipar
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 parametercj
.A Scilab function.
Its calling sequence must be
r=jac(t,y,ydot,cj)
and thejac
function must returnr=dg(t,y,ydot)/dy+cj*dg(t,y,ydot)/dydot
wherecj
is a real scalar.A list.
It must be as follows
list(jac,x1,x2,...)
where the calling sequence of the function
jac
is nowr=jac(t,y,ydot,cj,x1,x2,...)
jac
still returnsdg/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:
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 islist([],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 byjac
;r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)"
.Ifjac
returns a full matrix setinfo(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 setinfo(6)=0
.- info(7)
set
info(7)=1
ifydot0
is just an estimation,info(7)=0
ifg(t0,y0,ydot0)=0
.
- hd
a real vector which allows to store the
dassl
context and to resume integration.- r
a real matrix. Each column is the vector
[t;x(t);xdot(t)]
wheret
is time index for which the solution had been computed.
Description
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
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; 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');
See Also
Report an issue | ||
<< dasrt | Differential calculus, Integration | diff >> |