ode - solveur d'équations différentielles ordinaires
Ici on ne décrit l'usage de ode que pour des EDO explicites.
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).
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 .
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)] .
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" .
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.o',[],'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)
ode_discrete , ode_root , dassl , impl , odedc , odeoptions , csim , ltitr , rtitr ,