Please note that the recommended version of Scilab is 2024.1.0. This page might be outdated.
See the recommended documentation of this function
ode
solveur d'équations différentielles ordinaires
Séquence d'appel
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)
Paramètres
- y0
vecteur ou matrice réelle (conditions initiales).
- t0
réel (instant initial).
- t
vecteur réel (instants où la solution est renvoyée).
- f
fonction externe (fonction Scilab ou chaîne de caractères ou liste).
- type
une des chaînes de caractères :
"adams"
,"stiff"
,"rk"
,"rkf"
,"fix"
,"discrete"
,"roots"
- rtol, atol
constantes ou vecteurs réels de même taille que
y
.- jac
fonction externe (fonction Scilab, chaîne de caractères ou liste).
- w, iw
vecteurs réels
- ng
entier
- g
fonction externe (fonction Scilab, chaîne de caractères ou liste).
- k0
entier (instant initial).
- kvect
vecteur d'entiers.
Description
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 bibliothèques, 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 à :
- <aucun> :
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.- "adams" :
Problèmes non raides. Le solveur
lsode
du package ODEPACK est utilisé (schéma d'Adams).- "stiff" :
Pour les systèmes raides. Le solveur
lsode
du package ODEPACK est utilisé avec le schéma BDF.- "rk" :
Schéma de Runge-Kutta adaptatif d'ordre 4 (RK4).
- "rkf" :
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.
- "fix":
Identique à
"rkf"
, mais l'interface est simplifiée, i.e. uniquementrtol
etatol
sont communiqués au solveur.- "root":
Solveur d'EDO avec recherche de racines. Le solveur
lsodar
du package ODEPACK est utilisé. C'est une variante delsoda
permettant la recherche d'une racine d'une fonction vectorielle donnée. Voir ode_root pour plus de détails.- "discrete":
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, ett
est le vecteur des instants où l'on veut une approximation de la solution.y
est calculée ety
est la matricey=[y(t(1)),y(t(2)),...]
.Le paramètre d'entrée
f
deode
définit 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 procedurefex
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:
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 syntaxeydot = 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 matricep x q
au lieu d'un vecteur. Dans ce cas, on résout le système d'EDOn=p+q
dY/dt=F(t,Y)
oùY
est une matricep x q
. La condition initialeY0
doit aussi être une matricep x q
matrix et le résultat renvoyé parode
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
etatol
sont des valeurs seuil sur les erreurs estimées (relative et absolue) L'erreur estimée sury(i)
estrtol(i)*abs(y(i))+atol(i)
Si
rtol
et/ouatol
sont des constantesrtol(i)
et/ouatol(i)
prennent ces valeurs. Les valeurs par défaut dertol
etatol
sont respectivementrtol=1.d-5
etatol=1.d-7
pour la plupart des solveurs etrtol=1.d-3
etatol=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ètrejac
deode
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 êtreJ=jac(t,y)
où
t
est un scalaire (le temps) ety
un vecteur (l'état). La matriceJ
doit renvoyer df/dx i.e.J(k,i) = dfk /dxi
avecfk
= 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:
Dans la plupart des cas il n'est pas nécessaire d'utiliser
ml
,mu
etnrpd
, qui sont relatifs à la possibilité de stockage "bande" du JacobienSi
jac
est une liste, les mêmes conventions que pourf
s'appliquent.Les arguments optionnels
w
etiw
sont des vecteurs dans lesquels 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 appel précédent àode
.Plus d'options peuvent être passées aux solveurs d'ODEPACK en utilisant la variable
%ODEOPTIONS
. Voir odeoptions.
Exemples
// ---------- 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+'/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)
Voir aussi
- ode_discrete — ordinary differential equation solver, discrete time simulation
- ode_root — ordinary differential equation solver with root finding
- dassl — differential algebraic equation
- impl — differential algebraic equation
- odedc — discrete/continuous ode solver
- odeoptions — set options for ode solvers
- csim — simulation (time response) of linear system
- ltitr — discrete time response (state space)
- rtitr — discrete time response (transfer matrix)
- intg — intégration numérique adaptative
Bibliographie
Alan C. Hindmarsh, lsode and lsodi, two new initial value ordinary differential equation solvers, acm-signum newsletter, vol. 15, no. 4 (1980), pp. 10-11.
Fonctions Utilisées
Les sous-programmes associés se trouvent dans le répertoire SCI/modules/differential_equations/src/fortran/: lsode.f lsoda.f lsodar.f
Report an issue | ||
<< inttrap | Calcul différentiel, Intégration | bvode >> |