daskr
solveur de DAE avec traversées de zéros
Séquence d'appel
[r, nn [, hd]] = daskr(x0, t0, t [, atol [, rtol]], res [, jac], ng, surf [, info [, psol] [, pjac]] [, hd])
Paramètres
- x0
représente soit
y0
(ydot0
sera estimé pardaskr
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 deydot0
, assignezinfo(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 taille1
.atol
etrtol
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 dey
.- 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ésidur = g(t, y, ydot)
et un drapeau d'erreurires
.ires = 0
sires
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 maintenantr = res(t, y, ydot, x1, x2, ...)
res
retourne toujoursr = 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
etipar
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ètrecj
.Une fonction Scilab.
Sa séquence d'appel doit être
r = jac(t, y, ydot, cj)
et doit retournerr = 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ésormaisr = jac(t, y, ydot, cj, x1, x2,...)
jac
retourne toujoursdg/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 maintenantr = 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 estlist([], 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é parjac
;r(i - j + ml + mu + 1,j) = "dg(i)/dy(j)+cj*dg(i)/dydot(j)"
. Sijac
retourne une matrice pleine, fixerinfo(3)=[]
. Inutile siinfo(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 fixerinfo(6)=0
.- info(7)
si ydot0 est fixé tel que
g(t0, y0, ydot0) = 0
, alors fixerinfo(7)=[]
. Sinon, fixerinfo(7)=[+-1, ..., +-1]
, avecinfo(7)(i) = 1
si y(i) est une variable différentielle etinfo(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 routinepsol
si vous souhaitez que le solveur utilise des itérations de Krylov, sinon (méthode directe) fixerinfo(8)=0
.- info(9)
paramètres de Krylov. Inutile si vous avez fixé
info(8)=0
. Sinon, fixerinfo(9)=[]
ouinfo(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)=[]
. Fixerinfo(10)=1
si le solveur doit s'arrêter après le calcul des valeurs initiales, sinon fixerinfo(10)=0
.- info(11)
routine pour le calcul et la factorisation LU du préconditionneur pour
psol
. Inutile siinfo(8)=0
. Fixerinfo(11)=1
et fournir une routinepjac
si l'externalpsol
doit utiliser une routine spécifique, sinon fixerinfo(11)=0
.- info(12)
si vous souhaitez contrôler l'erreur localement sur toutes les variables, fixez
info(12)=[]
. Sinon, fixezinfo(12)=[+-1, ..., +-1]
, avecinfo(12)(i) = 1
si y(i) est une variable différentielle etinfo(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, fixerinfo(13)=[]
ouinfo(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
siinfo(8)=0
,15
sinon),- mxnj = nombre maximal d'évaluations du préconditionneur (par défaut
6
siinfo(8)=0
,2
sinon),- mxnh = nombre maximal de valeurs artificielles du pas
h
à tenter si info(7) ≠ [] (par défaut5
),- 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 fixerinfo(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é quepjac
a calculé auparavant et stocké danswp
etiwp
.Une fonction Scilab.
Sa séquence d'appel doit être
[r, ier] = psol(wp, iwp, b)
et doit retourner la solution du système dansr
et un drapeau d'erreurier
.Une liste.
Elle doit se présenter comme suit :
list(psol, x1, x2, ...)
où la séquence d'appel de
psol
estpsol(wp, iwp, b, x1, x2, ...)
psol
retourne toujours la solution dansr
.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 :
où les vecteursvoid 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)
wp
etiwp
contiennent le préconditionneur LU-factoriséP
,wp
représentant les valeurs etiwp
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ètrecj
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 vecteurswp
etiwp
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
estpjac(neq, t, y, ydot, h, cj, rewt, savr, x1, x2,...)
pjac
retourne toujoursdg/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
Historique
Version | Description |
5.5.0 | daskr ajouté |
Report an issue | ||
<< daeoptions | Intégration - dérivation | dasrt >> |