daskr

Paramètres

x0

représente soit y0 (ydot0 sera estimé par daskr avec zéro comme première estimation), soir la matrice [y0 ydot0]. g(t, y0, ydot0) doit être égal à zéro. Si vous ne connaissez qu'une estimation de ydot0, assignez info(7)=1.

y0

vecteur colonne réel des conditions initiales.

ydot0

vecteur colonne réel de la dérivée en temps de y à t0 (peut être une estimation.

t0

réel, temps initial.

t

réel, scalaire ou vecteur. Temps auxquels la solution est désirée. Notez que vous pouvez obtenir la solution à chaque étape de daskr en fixant info(2)=1.

atol, rtol

réels scalaires ou vecteurs colonnes de même taille que y ou tous deux de taille 1. atol et rtol représentent les tolérances d'erreur absolue et relative de la solution. Si ce sont des vecteurs, alors les tolérances sont spécifiées pour chaque composante de y.

res

external (fonction, liste ou chaîne de caractères). Calcule la valeur de g(t, y, ydot). Elle peut être :

• Une fonction Scilab.

  Sa séquence d'appel doit être [r, ires] = res(t, y, ydot) et doit retourner le résidu r = g(t, y, ydot) et un drapeau d'erreur ires. ires = 0 si res a correctement calculé r, ires = -1 si le résidu est localement non défini pour (t, y, ydot), ires = -2 si des paramètres sont hors du champ admissible.

• Une liste.

  Cette forme permet de passer des paramètres autres que t, y, ydot à la fonction. Elle doit se présenter comme suit :

  list(res, x1, x2, ...)

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

  r = res(t, y, ydot, x1, x2, ...)

  res retourne toujours r = g(t, y, ydot) comme fonction de (t, y, ydot, x1, x2, ...).

  Attention : cette forme ne doit pas être utilisée s'il n'y pas d'argument additionnel à passer à la fonction.

• Une chaîne de caractères.

  Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab.

  En C, la séquence d'appel doit être :

  void res(double *t, double y[], double yd[], double r[],
           int *ires, double rpar[], int ipar[])

  En Fortran, elle doit être :

  subroutine res(t, y, yd, r, ires, rpar, ipar)
  double precision t, y(*), yd(*),r(*),rpar(*)
  integer ires, ipar(*)

  Les tableaux rpar et ipar doivent être présents mais ne peuvent pas être utilisés.

jac

external (fonction, liste ou chaîne de caractères). Calcule la valeur de dg/dy + cj*dg/dydot pour une valeur donnée du paramètre cj.

• Une fonction Scilab.

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

• Une liste.

  Elle doit se présenter comme suit :

  list(jac, x1, x2, ...)

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

  r=jac(t, y, ydot, cj, x1, x2,...)

  jac retourne toujours dg/dy + cj*dg/dydot comme fonction de (t, y, ydot, cj, x1, x2,...).

• Une chaîne de caractères.

  Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab.

  En C, la séquence d'appel doit être :

  void jac(double *t, double y[], double yd[], double pd[],
           double *cj, double rpar[], int ipar[])

  En Fortran, elle doit être :

  subroutine jac(t, y, yd, pd, cj, rpar, ipar)
  double precision t, y(*), yd(*), pd(*), cj, rpar(*)
  integer ipar(*)

surf

external (fonction, liste ou chaîne de caractères). Calcule la valeur du vecteur colonne surf(t, y) à ng composantes. Chaque composante représente une surface. Elle doit être définie comme suit :

• Une fonction Scilab.

  Sa séquence d'appel doit être surf(t, y)

• Une liste.

  Elle doit se présenter comme suit :

  list(surf, x1, x2, ...)

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

  r = surf(t, y, x1, x2, ...)

• Une chaîne de caractères.

  Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab.

  En C, la séquence d'appel doit être :

  void surf(int *ny, double *t, double y[], int *ng, double gout[])

  En Fortran, elle doit être :

  subroutine surf(ny, t, y, ng, gout)
  double precision t, y(*), gout(*)
  integer ny, ng

info

liste contenant 14 éléments. La valeur par défaut est list([], 0, [], [], [], 0, [], 0, [], 0, 0, [], [], 1).

info(1)

réel scalaire donnant le temps maximal pour lequel g peut être évalué ou une matrice vide [] si aucune limite de temps n'est imposée.

info(2)

drapeau indiquant si daskr retourne ses valeurs intermédiaires calculées (= 1) ou seulement les temps indiqués par l'utilisateur (= 0).

info(3)

vecteur de deux éléments donnant la définition [ml,mu] de la matrice bande calculeé par jac; r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)". Si jac retourne une matrice pleine, fixer info(3)=[]. Inutile si info(8)=1.

info(4)

réel scalaire donnant la taille maximale du pas. Fixer info(4)=[] si illimité.

info(5)

réel scalaire donnant le pas initial. Fixer info(5)=[] si non spécifié.

info(6)

fixer info(6)=1 si la solution est non-négative, sinon fixer info(6)=0.

info(7)

si ydot0 est fixé tel que g(t0, y0, ydot0) = 0, alors fixer info(7)=[]. Sinon, fixer info(7)=[+-1, ..., +-1], avec info(7)(i) = 1 si y(i) est une variable différentielle et info(7)(i) = -1 si y(i) est une variable algébrique (si ses dérivées n'apparaissent pas explicitement dans la fonction g(t, y, ydot)).

info(8)

méthode directe / Krylov. Fixer info(8)=1 et founrnir une routine psol si vous souhaitez que le solveur utilise des itérations de Krylov, sinon (méthode directe) fixer info(8)=0.

info(9)

paramètres de Krylov. Inutile si vous avez fixé info(8)=0. Sinon, fixer info(9)=[] ou info(9)=[maxl kmp nrmax epli], où :

- maxl = nombre maximal d'itérations de l'algorithme GMRes (par défaut min(5, neq)),

- kmp = nombre de vecteurs sur lesquels l'orthogonalisation est faite dans GMRes (par défaut maxl),

- nrmax = nombre maximal de redémarrages de GMRes par intération non-linéaire (par défaut 5),

- epli = constante du test de convergence de GMRes (par défaut 0.05).

info(10)

conditions initiales. A ignorer si info(7)=[]. Fixer info(10)=1 si le solveur doit s'arrêter après le calcul des valeurs initiales, sinon fixer info(10)=0.

info(11)

routine pour le calcul et la factorisation LU du préconditionneur pour psol. Inutile si info(8)=0. Fixer info(11)=1 et fournir une routine pjac si l'external psol doit utiliser une routine spécifique, sinon fixer info(11)=0.

info(12)

si vous souhaitez contrôler l'erreur localement sur toutes les variables, fixez info(12)=[]. Sinon, fixez info(12)=[+-1, ..., +-1], avec info(12)(i) = 1 si y(i) est une variable différentielle et info(12)(i) = -1 si y(i) est une variable algébrique (si ses dérivées n'apparaissent pas explicitement dans la fonction g(t, y, ydot)).

info(13)

paramètres heuristiques. Ignorer si info(7)=[]. Sinon, fixer info(13)=[] ou info(13)=[mxnit mxnj mxnh lsoff stptol epinit], où :

- mxnit = nombre maximal d'itérations de Newton par évaluation du préconditionneur (par défaut 5 si info(8)=0, 15 sinon),

- mxnj = nombre maximal d'évaluations du préconditionneur (par défaut 6 si info(8)=0, 2 sinon),

- mxnh = nombre maximal de valeurs artificielles du pas h à tenter si info(7) ≠ [] (par défaut 5),

- lsoff = drapeau pour désactiver l'algorithme de recherche linéaire (lsoff = 0 pour activer, lsoff = 1 pour désactiver) (par défaut 0),

- stptol = pas minimal (dimmensionné) dans l'algorithme de recherche linéaire (par défaut (unit roundoff)^(2/3)),

- epinit = facteur déterminant dans le test de convergence de l'itération Newton (par défaut 0.01).

info(14)

verbosité. Fixer info(14)=1 pour une information minimale, info(14)=2 pour une information maximale, sinon fixer info(14)=0.

psol

external (fonction, liste ou chaîne de caractères). Résout un système linéraire P*x = b, où P est le préconditionneur LU-factorisé que pjac a calculé auparavant et stocké dans wp et iwp.

• Une fonction Scilab.

  Sa séquence d'appel doit être [r, ier] = psol(wp, iwp, b) et doit retourner la solution du système dans r et un drapeau d'erreur ier.

• Une liste.

  Elle doit se présenter comme suit :

  list(psol, x1, x2, ...)

  où la séquence d'appel de psol est

  psol(wp, iwp, b, x1, x2, ...)

  psol retourne toujours la solution dans r.

• Une chaîne de caractères.

  Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab

  En C, la séquence d'appel doit être :

  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)

  où les vecteurs wp et iwp contiennent le préconditionneur LU-factorisé P, wp représentant les valeurs et iwp les pivots utilisés dans la factorisation.

  En Fortran, elle doit être :

  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 (fonction, liste ou chaîne de caractères). Calcule la valeur de dg/dy + cj*dg/dydot pour une valeur donnée du paramètre cj et la LU-factorise en deux vecteurs, réel et entier.

• Une fonction Scilab.

  Sa séquence d'appel doit être [wp, iwp, ires] = pjac(neq, t, y, ydot, h, cj, rewt, savr) et en retour, les vecteurs wp et iwp doivent contenir toutes les informations du préconditionneur factorisé.

• Une liste.

  Elle doit se présenter comme suit :

  list(pjac, x1, x2, ...)

  où la séquence d'appel de pjac est

  pjac(neq, t, y, ydot, h, cj, rewt, savr, x1, x2,...)

  pjac retourne toujours dg/dy + cj*dg/dydot comme fonction de (neq, t, y, ydot, h, cj, rewt, savr, x1, x2, ...).

• Une chaîne de caractères.

  Elle doit se référer au nom d'une fonction C ou une routine Fortran reliée à Scilab

  En C, la séquence d'appel doit être :

  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)

  En Fortran, elle doit être :

  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

vecteur réel servant à stocker le contexte de daskr et reprendre l'intégration.

r

matrice réelle. Chaque colonne est le vecteur [t; x(t); xdot(t)] où t est l'indice en temps aulequel la solution a été calculée, x(t) est la valeur de la solution calculée, xdot(t) est la dérivée de la solution calculée.

nn

vecteur à deux entrées [times num], times est la valeur du temps auquel la surface est traversée, num est le nombre de surfaces traversées.

Description

Solution de l'équation différentielle implicite :

g(t, y, ydot) = 0
y(t0) = y0  et  ydot(t0) = ydot0

Retourne les temps de traversée de surface et le nombre de surfaces traversées dans nn.

Des exemples détaillés se trouvent dans SCI/modules/differential_equations/tests/unit_tests/daskr.tst

Exemples

// dy/dt = ((2*log(y)+8)/t-5)*y,  y(1) = 1,  1 <= t <= 6
// g1 = ((2*log(y)+8)/t-5)*y
// g2 = log(y) - 2.2491
y0 = 1; t = 2:6; t0 = 1; y0d = 3;
atol = 1.d-6; rtol = 0; ng = 2;

deff('[delta, ires] = res1(t, y, ydot)', 'ires = 0; delta = ydot-((2*log(y)+8)/t-5)*y')
deff('[rts] = gr1(t, y)', 'rts = [((2*log(y)+8)/t-5)*y; log(y)-2.2491]')

[yy, nn] = daskr([y0, y0d], t0, t, atol, rtol, res1, ng, gr1);
nn

// Retourne nn = [2.4698972 2]

Voir aussi

• ode — solveur d'équations différentielles ordinaires
• dasrt — DAE solver with zero crossing
• dassl — differential algebraic equation
• impl — differential algebraic equation
• fort — Fortran or C user routines call
• link — dynamic linker
• external — Objet Scilab, fonction externe ou routine

Historique