dae

Arguments

y0

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 .

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.

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.

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

psol

external (function, list or string). Solves a linear system P*x = b, with P being the factored preconditioner that routine pjac computed beforehand and stored in wp and iwp.

A Scilab function.

Its syntax must be [r, ier] = psol(wp, iwp, b) and must return the solution of the system in r and an error flag ier.
A list.

It must be as follows:
```
list(psol, x1, x2, ...)
```
where the syntax of psol is now
```
psol(wp, iwp, b, x1, x2, ...)
```
psol still returns the solution in r.

A character string.

It must refer to the name of a C function or a Fortran subroutine linked with Scilab

In C, the syntax must be:

void psol (int*neq, double*t, double*y, double*ydot, double*savr,
  double*wk, double*cj, double*wght, double*wp, int*iwp, double*b, double*eplin, int*ier, double*rpar, int*ipar)

where the arrays wp and iwp contain matrix elements of LU-factored preconditioner P, wp being the values and iwp the pivots used in the factorization.

In Fortran, it must be:

subroutine psol (neq, t, y, ydot, savr, wk, cj, wght,
  wp, iwp, b, eplin, ier, rpar, ipar)
double precision t,y(*), ydot(*), savr(*), wk(*), cj, wght(*), wp(*),
  b(*), eplin, rpar(*)
integer neq, iwp(*), ier, ipar(*)

pjac

external (function, list or string). Computes the value of dg/dy + cj*dg/dydot for a given value of parameter cj and LU-factorizes it in two arrays, real and integer.

A Scilab function.

Its syntax must be [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr) and in return, the arrays wp and iwp must contain all factored preconditioner information.
A list.

It must be as follows
```
list(pjac, x1, x2, ...)
```
where the syntax of pjac is
```
pjac(neq, t, y, ydot, h, cj, rewt, savr, x1, x2,...)
```
pjac still returns factorized dg/dy + cj*dg/dydot as a function of (neq, t, y, ydot, h, cj, rewt, savr, x1, x2, ...).

A character string.

It must refer to the name of a C function or a Fortran subroutine linked with Scilab

In C, the syntax must be:

void pjac (double*res, int*ires, int*neq, double*t, double*y, double*ydot, double*rewt, double*savr,
double*wk, double*h, double*cj, double*wp, int*iwp, int*ier, double*rpar, int*ipar)

In Fortran, it must be:

subroutine pjac (res, ires, neq, t, y, ydot, rewt, savr,
  wk, h, cj, wp, iwp, ier, rpar, ipar)
double precision res(*), t, y(*), ydot(*), rewt(*), savr(*),
  wk(*), h, cj, wp(*), rpar(*)
integer ires, neq, iwp(*), ier, ipar(*)

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 solves the differential algebraic equation and returns the evolution of y a given time points

g(t, x, xdot) = 0
x(t0) = x0 and xdot(t0) = xdot0

If xdot0 is not given in the y0 argument, the dae function tries to compute it solving g(t,x0,xdot0)=0.

If xdot0 is given in the y0 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.

The solvers

The type of problem solved and the method used depend on the value of the first optional argument type which can be one of the following strings:

<not given>:: the dassl solver.
root: the dasrt solver. See help on dae_root for more details.
root2:: the daskr solver. See help on dae_root for more details.

Examples

Example #1: use dassl solver

// 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: use dasrt solver with "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: use daskr solver with "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);

History

Версия	Описание
2024.1.0	dassl, dasrt and daskr become internal.

Report an issue
<< bvode	Дифференциальные уравнения	dae_root >>

Syntax

Arguments

Description

The solvers

Examples

See also

History