Please note that the recommended version of Scilab is 2024.1.0. This page might be outdated.
See the recommended documentation of this function
dae
Differential algebraic equations solver
Syntax
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]) [y, rd] = dae("root2", initial, t0, t, res, ng, surface) [y, rd [,hd]] = dae("root2", initial, t0, t [[,rtol], atol], res [,jac], ng, surface [, psol, pjac] [, hd])
Arguments
- initial
a column vector. It may be equal to
x0
or[x0;xdot0]
. Wherex0
is the state value at initial timet0
andxdot0
is the initial state derivative value or an estimation of it (see below).- t0
a real number, the initial time.
- t
a 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
.- atol
a real scalar or a column vector of same size as
x0
, the absolute error tolerance of solution. Ifatol
is a vector the tolerances are specified for each component of the state.- rtol
a real scalar or a column vector of same size as
x0
, the relative error tolerance of solution. Ifrtol
is a vector the tolerances are specified for each component of the state.- res
an external function, computes the value of
g(t,y,ydot)
. It may be- a Scilab function
In this case, the syntax 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
.ires = -1
if residue is locally not defined forg(t,x,xdot)
.ires =-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:
list(res,p1,p2,...)
where the syntax of the function
res
is nowr=res(t,y,ydot,p1,p2,...)
res
still returns the residual value as a function of(t,x,xdot,x1,x2,...)
, andp1, 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,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 syntax 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 valuex
the state arrayxdot
the array of state derivativesres
the array of residualsires
the execution indicatorrpar
is the array of floating point parameter values, needed but cannot be set by thedae
functionipar
is the array of floating integer parameter values, needed but cannot be set by thedae
function
- jac
an external, computes the value of
dg/dx+cj*dg/dxdot
for a given value of parametercj
. It may be- a Scilab function
Its syntax must be
r=jac(t,x,xdot,cj)
and thejac
function must returnr=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot
wherecj
is a real scalar.- a list
This form of external is used to pass parameters to the function. It must be as follows:
list(jac,p1,p2,...)
where the syntax of the function
jac
is nowr=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,
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 syntax 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
- surface
an external, computes the value of the column vector
surface(t,x)
withng
components. Each component defines a surface.- a Scilab function
Its syntax must be
r=surface(t,x)
, this function must return a vector withng
elements.- a list
This form of external is used to pass parameters to the function. It must be as follows:
list(surface,p1,p2,...)
where the syntax of the function
surface
is nowr=surface(t,x,p1,p2,...)
- a character string
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
<s_name>(nx, t, x, ng, r, rpar, ipar)
double precision t, x(*), r(*), rpar(*)
integer nx, ng,ipar(*)
The C syntax must be
C2F(<s_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 andr
is the results array.
- rd
a vector with two entries
[times num]
wheretimes
is the value of the time at which the surface is crossed,num
is the number of the crossed surface- pjac
external (function, list or string). See daskr
- psol
external (function, list or string). See daskr
- hd
a real vector, as an output it stores the
dae
context. It can be used as an input argument to resume integration (hot restart).- y
a real matrix. If
%DAEOPTIONS(2)=1
, each column is the vector[t;x(t);xdot(t)]
wheret
is time index for which the solution has been computed. Elsey
is the vector[x(t);xdot(t)]
.
Description
The dae
function is a gateway built above the
dassl, dasrt
and daskr
functions designed for implicit differential equations integration.
Option "root"
calls the dasrt routine,
and "root2"
calls daskr.
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
argument 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
and
modules/differential_equations/tests/unit_tests/daskr.tst
Examples
Example #1: dassl (no root finding)
// 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 #2: dasrt ("root")
// Example with C code (C compiler needed) //-------------------------------------------------- bOK = haveacompiler(); if bOK <> %t [btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info'); return end //-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];}'] previous_dir = pwd(); cd TMPDIR; mputl(code, 't22.c') //-2- Compile and load them ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce'); exec('t22loader.sce') //-3- Run rtol = [1.d-6; 1.d-6]; atol = [1.d-6; 1.d-4]; t0 = 0; t = [20:20:200]; y0 = [2; 0]; y0d = [0; -2]; 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))); cd(previous_dir);
Example #3: daskr ("root2"), using default 'psol' and 'pjac' routines
// Example with C code (C compiler needed) //-------------------------------------------------- bOK = haveacompiler(); if bOK <> %t [btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info'); return end //-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];}'] previous_dir = pwd(); cd TMPDIR; mputl(code, 't22.c') //-2- Compile and load them ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce'); exec('t22loader.sce') //-3- Run rtol = [1.d-6; 1.d-6]; atol = [1.d-6; 1.d-4]; t0 = 0; t = [20:20:200]; y0 = [2; 0]; y0d = [0; -2]; 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 %DAEOPTIONS = list([] , 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1); [yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1'); 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("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd); plot(yy(1, 1), yy(2, 1), 'r+') xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1))); cd(previous_dir);
See also
- ode — ordinary differential equation solver
- daeoptions — set options for dae solver
- dassl — differential algebraic equation
- dasrt — DAE solver with zero crossing
- daskr — DAE solver with zero crossing
- impl — differential algebraic equation
- call — Fortran or C user routines call
- link — dynamic linker
- external — Scilab Object, external function or routine
Report an issue | ||
<< bvode | Differential calculus, Integration | daeoptions >> |