dae

Paramètres

initial

vecteur colonne. Peut valoir x0 ou [x0;xdot0], où x0 est la valeur de l'état au temps initial t0 et xdot0 est la valeur (ou une estimation) de l'état dérivé au temps initial (voir ci-après).

t0

réel, le temps initial.

t

scalaire ou vecteur réel. Les instants auxquels la solution est voulue. La solution peut s'obtenir à chaque étape de dae en initialisant %DAEOPTIONS(2)=1 .

atol

un scalaire ou vecteur colonne réel de même taille que x0, la tolérance absolue de la solution. Si atol est un vecteur, les tolérances sont spécifiées pour chaque composante de l'état.

rtol

un scalaire ou vecteur colonne réel de même taille que x0, la tolérance relative de la solution. Si atol est un vecteur, les tolérances sont spécifiées pour chaque composante de l'état.

res

une fonction external, calcule la valeur de g(t,y,ydot). Peut être

une fonction Scilab

Dans ce cas, sa séquence d'appel doit être [r,ires]=res(t,x,xdot) et res doit retourner le résidu r=g(t,x,xdot) et un drapeau d'erreur ires.

ires = 0 si res est bien parvenu à calculer r.

ires = -1 si le résidu est localement non-défini pour g(t,x,xdot).

ires =-2 si le paramètres ne sont pas admissibles.

une liste

Cette forme d'external sert à passer des paramètres à la fonction. Elle doit se présenter comme suit:

list(res,p1,p2,...)

où la séquence d'appel de la fonction res est

r=res(t,y,ydot,p1,p2,...)

res retourne toujours le résidu comme fonction de (t,x,xdot,x1,x2,...), et p1, p2,... sont des paramètres de la fonction.

une chaîne de caractères

Elle doit se référer au nom d'une routine C ou Fortran. Supposant que <r_name> est le nom donné,

• La séquence d'appel en Fortran doit être

  <r_name>(t,x,xdot,res,ires,rpar,ipar)

  double precision t,x(*),xdot(*),res(*),rpar(*)

  integer ires,ipar(*)

• La séquence d'appel en C doit être

  C2F(<r_name>)(double *t, double *x, double *xdot, double *res, integer *ires, double *rpar, integer *ipar)

où

• t est la valeur actuelle du temps

• x est la valeur de l'état

• xdot est la valeur de l'état dérivé

• res la valeur du résidu

• ires indicateur de complétion

• rpar est un vecteur de paramètres nécéssaires mais non initialisables par dae.

• ipar est un vecteur de paramètres entiers nécéssaires mais non initialisables par dae

jac

une fonction external, calcule la valeur de dg/dx+cj*dg/dxdot pour un paramètre cj donné. Elle peut être:

une fonction Scilab

Sa séquence d'appel doit être r=jac(t,x,xdot,cj) et la fonction jac doit retourner r=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot où cj est un scalaire réel.

une liste

Cette forme d'external est utilisée pour passer des paramètres à la fonction. Elle doit se présenter comme suit :

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

où la séquence d'appel de la fonction jac est

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

jac retourne toujours dg/dx+cj*dg/dxdot comme fonction de (t,x,xdot,cj,p1,p2,...).

une chaîne de caractères

elle doit se référer au nom d'une routine C ou Fortran. Supposant que <j_name> est le nom donné,

• La séquence d'appel en Fortran doit être

  <j_name>(t, x, xdot, r, cj, ires, rpar, ipar)

  double precision t, x(*), xdot(*), r(*), ci, rpar(*)

  integer ires, ipar(*)

• La séquence d'appel en C doit être

  C2F(<j_name>)(double *t, double *x, double *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)

où t, x, xdot, ires, rpar, ipar ont la même définition qu'au-dessus, r est la valeur du résultat

surface

une fonction external, calcule a valeur du vecteur colonne surface(t,x) à ng éléments. Chaque élément définit une surface. Elle peut être

une fonction Scilab

Sa séquence d'appel doit être r=surface(t,x), cette fonction doit retourner un vecteur colonne à ng éléments.

une liste

Cette forme d'external est utilisé pour passer des paramètres à la fonction. Elle doit se présenter comme suit :

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

où la qéquence d'appel de la fonction surface est

r=surface(t,x,p1,p2,...)

une chaîne de caractères

elle doit se référer au nom d'une routine C ou Fortran. Supposant que <s_name> est le nom donné,

• La séquence d'appel en Fortran doit être

  <s_name>(nx, t, x, ng, r, rpar, ipar)

  double precision t, x(*), r(*), rpar(*)

  integer nx, ng,ipar(*)

• La séquence d'appel en C doit être

  C2F(<s_name>)(double *t, double *x, double *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)

où t, x, rpar, ipar ont les mêmes définitions qu'au-dessus, ng est le nombre de surfaces, nx est la dimension de l'état et r est le vecteur résultat.

rd

un vecteur à deux entrées [times num] où times est la valeur des temps auquel les surfaces sont traversées, num est le nombre de surfaces traversées.

pjac

external (fonction, liste or chaîne de caractères). Voir daskr

psol

external (fonction, liste or chaîne de caractères). Voir daskr

hd

vecteur réel, stocke en sortie le contexte de dae. Peut être utilisé comme paramètre d'entrée pour reprendre l'intégration (reprise à chaud).

y

matrice réelle. Si %DAEOPTIONS(2)=1 , chaque colonne est le vecteur [t;x(t);xdot(t)] où t est le temps auquel la solution à été calculée. Sinon y est le vecteur [x(t);xdot(t)].

Description

La fonction dae est une gateway construite sur les solveurs dassl, dasrt et daskr destinés à l'intégration implicite d'équations différentielles.

L'option "root" appelle la routine dasrt, et "root2" appelle daskr.

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

Si xdot0 n'est pas donné en paramètre initial, dae tente de la calculer en résolvant g(t,x0,xdot0)=0.

Si xdot0 est donné en paramètre initial, il peut soit satifaire g(t,x0,xdot0)=0, soit en être une approximation. Dans le dernier cas, %DAEOPTIONS(7) doit être initialisé à 1.

Des exemples détaillés utilisant des externals codées en C et en Scilab se trouvent dans modules/differential_equations/tests/unit_tests/dassldasrt.tst et modules/differential_equations/tests/unit_tests/daskr.tst

Exemples

Exemple #1: dassl (pas de traversée de surface)

// Exemple avec du code Scilab
//--------------------------------------------------
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); // Retourne les points d'observation demandés

%DAEOPTIONS = list([], 1, [], [], [], 0, 0); // Demande à dae les points à retourner
y = dae([x0, xd0], 0, 4d10, chemres); // Sans jacobian
y = dae([x0, xd0], 0, 4d10, chemres, chemjac); // Avec jacobien

Exemple #2: dasrt ("root")

// Exemple avec du code C (compilateur C requis)
//--------------------------------------------------
bOK = haveacompiler();
if bOK <> %t
    [btn] = messagebox(["Vous avez besoin d''un compilateur C pour cet exemple."; "Arrêt de l''exécution."], "Problème de Software", "info");
    return
end

//-1- Crée les codes C dans TMPDIR - équation de Vanderpol, forme implicite
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 et charge
ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
exec('t22loader.sce')

//-3- Exécute
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;

// Simulation simple
t = 0:0.003:300;
yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
clf(); plot(yy(1, :), yy(2, :))
// Trouve le premier point où 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)));

// Reprise à chaud pour le point suivant
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);

Exemple #3: daskr ("root2"), avec les routines 'psol' et 'pjac' par défaut

// Exemple avec du code C (compilateur C requis)
//--------------------------------------------------
bOK = haveacompiler();
if bOK <> %t
    [btn] = messagebox(["Vous avez besoin d''un compilateur C pour cet exemple."; "Arrêt de l''exécution."], "Problème de Software", "info");
    return
end

//-1- Crée les codes C dans TMPDIR - équation de Vanderpol, forme implicite
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 et charge
ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
exec('t22loader.sce')

//-3- Exécute
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;

// Simulation simple
t = 0:0.003:300;
yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
clf(); plot(yy(1, :), yy(2, :))
// Trouve le premier point où 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)));

// Reprise à chaud pour le point suivant
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);

Voir aussi

• ode — solveur d'équations différentielles ordinaires
• daeoptions — set options for dae solver
• dassl — differential algebraic equation
• dasrt — DAE solver with zero crossing
• daskr — solveur de DAE avec traversées de zéros
• impl — differential algebraic equation
• call — Fortran or C user routines call
• link — dynamic linker
• external — Objet Scilab, fonction externe ou routine