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+qdY/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,