dae - Differential algebraic equations solver (original) (raw)

Scilab 5.3.3

Please note that the recommended version of Scilab is 2026.0.1. This page might be outdated.
See the recommended documentation of this function

Scilab help >> Differential Equations, Integration > 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 x0or[x0;xdot0]. Where x0 is the state value at initial time t0 andydot0 is the initial state derivative value or an estimation of it (see below).

t0

a real number, the initial time.

t

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](daeoptions.html)(2)=1.

rtol

a real scalar or a column vector of same size asx0. 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 asx0. The absolute error tolerance of solution. If atol is a vector the tolerances are specified for each component of the state.

res

an external. Computes the value ofg(t,y,ydot). It may be

a Scilab function

In this case, Its calling sequence must be[r,ires]=res(t,x,xdot) andres must return the residuer=g(t,x,xdot) and error flagires. ires = 0 ifres succeeds to computer, =-1 if residue is locally not defined for (t,x,xdot),=-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:

where the calling sequence of the functionres 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.

where

jac

an external. Computes the value ofdg/dx+cj*dg/dxdot for a given value of parametercj. It may be

a Scilab function

Its calling sequence must ber=jac(t,x,xdot,cj) and thejac function must returnr=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:

where the calling sequence of the functionjac is now

r=jac(t,x,xdot,p1,p2,...)

jac still returnsdg/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,

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 vectorsurface(t,x) with ng components. Each component defines a surface.

a Scilab function

Its calling sequence must ber=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:

where the calling sequence of the functionsurface is now

character string

it must refer to the name of a C or fortran routine. Assuming that <s_name> is the given name,

where t, x, rpar, ipar have similar definition as above, ngis the number of surfaces, nx the dimension of the state and r is the results array.

rd

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

hd

a real vector, an an output it stores thedae context. It can be used as an input argument to resume integration (hot restart).

y

real matrix . If [%DAEOPTIONS](daeoptions.html)(2)=1, each column is the vector [t;x(t);xdot(t)] wheret 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 thedassl 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 theinitial argument, the dae function tries to compute it solving g(t,x0,xdot0)=0,

if xdot0 is given in theinitial argumente it may be either a compatible derivative satisfying g(t,x0,xdot0)=0 or an approximate value . In the latter case[%DAEOPTIONS](daeoptions.html)(7) must be set to 1.

Detailed examples using Scilab and C coded externals are given inmodules/differential_equations/tests/unit_tests/dassldasrt.tst

Examples

function [r, ires]=chemres(t, y, yd) r(1) = -0.04y(1) + 1d4y(2)y(3) - yd(1); r(2) = 0.04y(1) - 1d4y(2)y(3) - 3d7y(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 , 1d4y(3) , 1d4y(2); 0.04 ,-1d4y(3)-23d7y(2)-cj ,-1d4y(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);

%DAEOPTIONS=list([],1,[],[],[],0,0); y=dae([x0,xd0],0,4d10,chemres); y=dae([x0,xd0],0,4d10,chemres,chemjac);

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.0y[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')

ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/t22loader.sce'); exec(TMPDIR+'/t22loader.sce')

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;

t=0:0.003:300; yy=dae([y0,y0d],t0,t,atol,rtol,'res22','jac22'); clf();plot(yy(1,:),yy(2,:))

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

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

See Also