dae

Arguments

initial

a column vector. It may be equal to x0 or [x0;xdot0]. Where x0 is the state value at initial time t0 and xdot0 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. If atol 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. If rtol 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) and res must return the residue r=g(t,x,xdot) and error flag ires.

ires = 0 if res succeeds to compute r.

ires = -1 if residue is locally not defined for g(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 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,

• 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 value

• x the state array

• xdot the array of state derivatives

• res the array of residuals

• ires the execution indicator

• rpar is the array of floating point parameter values, needed but cannot be set by the dae function

• ipar is the array of floating integer parameter values, needed but cannot be set by the dae function

jac

an external, computes the value of dg/dx+cj*dg/dxdot for a given value of parameter cj. It may be

a Scilab function

Its syntax must be r=jac(t,x,xdot,cj) and the jac function must return r=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:

list(jac,p1,p2,...)

where the syntax of the function jac is now

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

jac still returns dg/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) with ng components. Each component defines a surface.

a Scilab function

Its syntax must be r=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:

list(surface,p1,p2,...)

where the syntax of the function surface is now

r = 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 and r is the results array.

rd

a vector with two entries [times num] where times 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)] where t is time index for which the solution has been computed. Else y 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