Scilab Function

ode - solveur d'équations différentielles ordinaires

### Sequence 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)`

### Parametres

- `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 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 à :

- <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. uniquement `rtol` et `atol` 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 de `lsoda` 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, 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 `f` 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 `f` est une fonction Scilab, sa syntaxe doit être :

    ydot = f(t,y)

où `t` est un scalaire (le temps) et `y` un vecteur (l'état). Cette fonction renvoie le second membre de l'équation différentielle dy/dt=f(t,y).

Si `f` est une chaîne de caractères, elle désigne le nom d'une subroutine Fortran ou C, i.e. si `ode(y0,t0,t,"fex")` est la commande, alors la subroutine `fex` est appelée. Cette routine doit avoir la liste d'appel suivante : `f(n,t,y,ydot)`. La routine doit être liée dynamiquement à Scilab avec la fonction `link`. Voir les exemples dans les fichiers `SCIDIR/routines/default/README` et `SCIDIR/routines/default/Ex-ode.f`.

L'argument `f` peut aussi être une liste : si `ode(y0,t0,t,lst)` est la commande, alors `lst` doit être une liste avec la structure suivante :

    lst=list(f,u1,u2,...un)

où `f` est une fonction avec la syntaxe :

    ydot = f(t,y,u1,u2,...,un)

cela permet de passer des paramètres sous forme d'arguments supplémentaires de `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. Cette routine doit avoir la liste d'appel suivante : `jac(n,t,y,ml,mu,J,nrpd)`. Dans la plupart des cas il n'est pas nécessaire d'utiliser `ml`, `mu` et `nrpd` (voir les exemples dans `SCIDIR/routines/default/Ex-ode.f`).

Si `jac` est une liste, les mêmes conventions que pour `f` s'appliquent.

Les arguments optionnels `w` et `iw` sont des vecteurs permettant de redémarrer l'intégration au point où elle s'était arrêtée à la sortie de ode.

Plus d'options peuvent être passées aux solveurs d'ODEPACK en utilisant la variable `%ODEOPTIONS`. Voir le help de odeoptions.

### Exemples

    // EDO à une dimension
    // dy/dt=y^2-y sin(t)+cos(t), y(0)=0
    deff("[ydot]=f(t,y)","ydot=y^2-y*sin(t)+cos(t)")
    y0=0;t0=0;t=0:0.1:%pi;
    y=ode(y0,t0,t,f)
    plot(t,y)
    // Simulation de dx/dt = A x(t) + B u(t) avec u(t)=sin(omega*t),
    // x0=[1;0]
    // la solution x(t) est désirée en t=0.1, 0.2, 0.5 ,1.
    // A et u sont passées dans une liste
    // et B et omega sont des variables globales
    deff("[xdot]=linear(t,x,A,u)","xdot=A*x+B*u(t)")
    deff("[ut]=u(t)","ut=sin(omega*t)")
    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))
    //
    // EDO matricielle
    // Equation différentielle de Ricatti
    // Xdot=A'*X + X*A - X'*B*X + C , X(0)=Identité
    // Solution en t=[1,2]
    deff("[Xdot]=ric(t,X)","Xdot=A''*X+X*A-X''*B*X+C")
    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)
    A=[1,1;0,2];
    deff("[xdot]=f(t,x)","xdot=A*x");
    ode(eye(A),0,1,f)
    ode("adams",eye(A),0,1,f)
    // EDO raide, avec la jacobienne fournie
    A=[10,0;0,-1];
    deff("[xdot]=f(t,x)","xdot=A*x");
    deff("[J]=Jacobian(t,y)","J=A")
    ode("stiff",[0;1],0,1,f,Jacobian)

### Voir aussi

`ode_discrete` `ode_root` `dassl` `impl` `odedc` `odeoptions` `csim` `ltitr` `rtitr`