Scilab Website | Contribute with GitLab | Mailing list archives | ATOMS toolboxes
Scilab Online Help
2024.1.0 - 日本語


dae

Differential algebraic equations solver

Syntax

y = dae(y0, t0, t, res)
[y [,hd]] = dae(y0, t0, t [[,rtol], atol], res [,jac] [,hd])
[y, rd] = dae("root", y0, t0, t, res, ng, surface)
[y, rd [,hd]] = dae("root", y0, t0, t [[,rtol], atol], res [,jac], ng, surface [,hd])
[y, rd] = dae("root2", y0, t0, t, res, ng, surface)
[y, rd [,hd]] = dae("root2", y0, t0, t [[,rtol], atol], res [,jac], ng, surface [, psol, pjac] [, hd])

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

See also

  • dae_root — differential algebraic equation solver with roots finding
  • ode — 常微分方程式ソルバ
  • daeoptions — daeソルバのオプションを設定
  • impl — 微分代数方程式
  • call — Fortran or C user routines call
  • link — dynamic linker
  • external — Scilabオブジェクト, 外部関数またはルーチン

History

バージョン記述
2024.1.0 dassl, dasrt and daskr become internal.
Report an issue
<< bvode 微分方程式 dae_root >>
Copyright (c) 2022-2024 (Dassault Systèmes)
Copyright (c) 2017-2022 (ESI Group)
Copyright (c) 2011-2017 (Scilab Enterprises)
Copyright (c) 1989-2012 (INRIA)
Copyright (c) 1989-2007 (ENPC)
with contributors 		Last updated:
Mon Jun 17 17:54:17 CEST 2024