ode — solveur d'équations différentielles ordinaires
y=ode(y0,t0,t,f)
[y,w,iw]=ode([type],y0,t0,t [,rtol [,atol]],f [,jac] [,w,iw])
[y,rd,w,iw]=ode("root",y0,t0,t [,rtol [,atol]],f [,jac],ng,g [,w,iw])
y=ode("discrete",y0,k0,kvect,f)vecteur ou matrice réelle (conditions initiales).
réel (instant initial).
vecteur réel (instants où la solution est renvoyée).
fonction externe (fonction Scilab ou chaîne de caractères ou liste).
une des chaînes de caractères : "adams" "stiff" "rk"
"rkf" "fix" "discrete" "roots"
constantes ou vecteurs réels de même taille que
y.
fonction externe (fonction Scilab, chaîne de caractères ou liste).
vecteurs réels
entier
fonction externe (fonction Scilab, chaîne de caractères ou liste).
entier (instant initial).
vecteur d'entiers.
ode est la fonction utilisée pour approcher la
solution d'une équation différentielle ordinaire (EDO) explicite du
premier ordre en temps, définie par : dy/dt=f(t,y) , y(t0)=y0. Il s'agit
d'une interface vers diverses librairies, en particulier ODEPACK. Le type
du problème et la méthode utilisée dépendent de la valeur du premier
argument optionnel type qui peut être égal à :
le solveur lsoda du package ODEPACK est
utilisé par défaut. Il choisit automatiquement entre un schéma
prédicteur-correcteur d'Adams et un schéma adapté au systèmes raides
(stiff) de type "Backward Differentiation Formula" (BDF).
Initialement le schéma adapté aux système non raides est choisi puis
la méthode adaptée est ensuite choisie dynamiquement.
Problèmes non raides. Le solveur lsode du
package ODEPACK est utilisé (schéma d'Adams).
Pour les systèmes raides. Le solveur lsode
du package ODEPACK est utilisé avec le schéma BDF.
Schéma de Runge-Kutta adaptatif d'ordre 4 (RK4).
Formules de Shampine et Watts basées sur les paires de Runge-Kutta Fehlberg d'ordre 4 et 5 (RKF45). Bien pour les problèmes non raides ou moyennement raides, lorsque le calcul du second membre n'est pas trop coûteux. Cette méthode est à éviter si l'on recherche une très grande précision.
Identique à "rkf", mais l'interface est simplifiée, i.e.
uniquement rtol et atol sont
communiqués au solveur.
Solveur d'EDO avec recherche de racines. Le solveur
lsodar du package ODEPACK est utilisé. C'est une
variante de lsoda permettant la recherche d'une
racine d'une fonction vectorielle donnée. Voir ode_root pour plus de
détails.
Simulation en temps discret. Voir ode_discrete pour plus de détails.
Ici on ne décrit l'usage de ode que pour des EDO
explicites.
L'appel le plus simple de ode est du type :
y=ode(y0,t0,t,f) où y0 est le
vecteur des conditions initiales, t0 est le temps
initial, et t est le vecteur des instants où l'on
veut une approximation de la solution. y est
calculée et y est la matrice
y=[y(t(1)),y(t(2)),...].
Le paramètre d'entrée f de
ode défini le membre de droite de léquation
différentielle du premier ordre dy/dt=f(t,y). C'est un external qui
peut être :
Soit une fonction Scilab, sa syntaxe doit être ydot
= f(t,y) où t est un scalaire (le
temps), y un vecteur (l'état). Cette fonction
renvoie le second membre de l'équation différentielle
dy/dt=f(t,y).
Soit une chaîne de caractères, elle désigne le nom d'une
subroutine Fortran ou une procédure C, i.e. si
ode(y0,t0,t,"fex") est la commande, alors la
procedure fex est appelée.
Si c'est une subroutine Fortran, sa liste d'appel doit être
subroutine fex(n,t,y,ydot) integer n double precision t,y(*),ydot(*)
Si c'est une fonction C son prototype doit être:
void fex(int *n,double *t,double *y,double *ydot)
Cet external peut être compilé par l'utilitaire ilib_for_link et chargé dynamiquement par la fonction link.
Soit une liste avec la structure suivante
list(vrai_f,u1,u2,...un) où
vrai_f est une fonction avec la syntaxe
ydot = vrai_f(t,y,u1,u2,...,un)
Cette syntaxe permet de passer des paramètres sous forme
d'arguments supplémentaires de vrai_f.
La fonction f peut renvoyer une matrice
p x q au lieu d'un vecteur. Dans ce cas, on résout
le système d'EDO n=p+q
dY/dt=F(t,Y) où Y est une
matrice p x q. La condition initiale
Y0 doit aussi être une matrice p x
q matrix et le résultat renvoyé par ode
est la matrice: p x q(T+1) égale à
[Y(t_0),Y(t_1),...,Y(t_T)].
Des paramètres optionnels contrôlent la tolérance du schéma :
rtol et atol sont des valeurs
seuil sur les erreurs estimées (relative et absolue) L'erreur estimée
sur y(i) est
rtol(i)*abs(y(i))+atol(i)
Si rtol et/ou atol sont
des constantes rtol(i) et/ou
atol(i) prennent ces valeurs. Les valeurs par
défaut de rtol et atol sont
respectivement rtol=1.d-5 et
atol=1.d-7 pour la plupart des solveurs et
rtol=1.d-3 et atol=1.d-4 pour
"rfk" et "fix".
Pour les problèmes raides, il est recommandé de fournir la
jacobienne du second membre sous forme de l'argument optionnel
jac. Le paramètre jac de
ode est par exemple une fonction Scilab, dont la
syntaxe est imposée, ou le nom d'une subroutine Fortran ou C (chaîne
de caractères) ou une liste.
Si jac est une fonction Scilab sa syntaxe
doit être J=jac(t,y)
où t est un scalaire (le temps) et
y un vecteur (l'état). La matrice
J doit renvoyer df/dx i.e. J(k,i) = dfk
/dxi avec fk = k-ième composante de
f.
Si f est une chaîne de caractères, elle
désigne le nom d'une subroutine Fortran ou C.
En Fortran, Cette routine doit avoir la liste d'appel suivante :
subroutine fex(n,t,y,ml,mu,J,nrpd) integer n,ml,mu,nrpd double precision t,y(*),J(*)
Si c'est une fonction C son prototype doit être:
void fex(int *n,double *t,double *y,int *ml,int *mu,double *J,int *nrpd,)
Dans la plupart des cas il n'est pas nécessaire d'utiliser
ml, mu et
nrpd, qui sont relatif aà la possibilité de
stockage "bande" du Jacobien
Si jac est une liste, les mêmes conventions
que pour f s'appliquent.
Les arguments optionnels w et
iw sont des vecteurs ou le solveur stocke des
informations sur son état(voir ode_optional_output pour plus de
détails) . Lorsque ces paramêtres sont utilisés comme argument
d'entrée, ils permettent de redémarrer l'intégration au point où elle
s'était arrêtée à la sortie d'un apple précédent à
ode.
Plus d'options peuvent être passées aux solveurs d'ODEPACK en
utilisant la variable %ODEOPTIONS. Voir odeoptions.
// ---------- EDO Simple (external : fonction Scilab)
// dy/dt=y^2-y sin(t)+cos(t), y(0)=0
function ydot=f(t,y),ydot=y^2-y*sin(t)+cos(t),endfunction
y0=0;t0=0;t=0:0.1:%pi;
y=ode(y0,t0,t,f)
plot(t,y)
// ---------- EDO Simple (external : code C)
ccode=['#include <math.h>'
'void myode(int *n,double *t,double *y,double *ydot)'
'{'
' ydot[0]=y[0]*y[0]-y[0]*sin(*t)+cos(*t);'
'}']
mputl(ccode,TMPDIR+'/myode.c') //create the C file
ilib_for_link('myode','myode.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/loader.sce');//compile
exec(TMPDIR+'/loader.sce') //incremental linking
y0=0;t0=0;t=0:0.1:%pi;
y=ode(y0,t0,t,'myode');
// ---------- Simulation de dx/dt = A x(t) + B u(t) avec u(t)=sin(omega*t),
// x0=[1;0]
// solution x(t) desired at t=0.1, 0.2, 0.5 ,1.
// A and u function are passed to RHS function in a list.
// B and omega are passed as global variables
function xdot=linear(t,x,A,u),xdot=A*x+B*u(t),endfunction
function ut=u(t),ut=sin(omega*t),endfunction
A=[1 1;0 2];B=[1;1];omega=5;
ode([1;0],0,[0.1,0.2,0.5,1],list(linear,A,u))
// ----------Integration de l'équation différentielle de Riccati (état matriciel)
// Xdot=A'*X + X*A - X'*B*X + C , X(0)=Identity
// Solution at t=[1,2]
function Xdot=ric(t,X),Xdot=A'*X+X*A-X'*B*X+C,endfunction
A=[1,1;0,2]; B=[1,0;0,1]; C=[1,0;0,1];
t0=0;t=0:0.1:%pi;
X=ode(eye(A),0,t,ric)
// ---------- Calcul de exp(A) (état matriciel)
A=[1,1;0,2];
function xdot=f(t,x),xdot=A*x;,endfunction
ode(eye(A),0,1,f)
ode("adams",eye(A),0,1,f)
// ---------- Calcul de exp(A) (état matriciel, cas raide, jacobien fourni)
A=[10,0;0,-1];
function xdot=f(t,x),xdot=A*x,endfunction
function J=Jacobian(t,y),J=A,endfunction
ode("stiff",[0;1],0,1,f,Jacobian)