Differential_equations help: fix dae tolerances
[scilab.git] / scilab / modules / differential_equations / help / fr_FR / dae.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) 2013 - Scilab Enterprises - Paul Bignier : added "roots2" (daskr)
5  * Copyright (C) 2008 - INRIA
6  * ...
7  *
8  * This file must be used under the terms of the CeCILL.
9  * This source file is licensed as described in the file COPYING, which
10  * you should have received as part of this distribution.  The terms
11  * are also available at
12  * http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
13  *
14  -->
15 <refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:ns3="http://www.w3.org/1999/xhtml" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" xmlns:scilab="http://www.scilab.org" xml:id="dae" xml:lang="en">
16     <refnamediv>
17         <refname>dae</refname>
18         <refpurpose>Solveur de système d'Equations Différentielles Algébriques</refpurpose>
19     </refnamediv>
20     <refsynopsisdiv>
21         <title>Séquence d'appel</title>
22         <synopsis> y = dae(initial, t0, t, res)
23             [y [,hd]] = dae(initial, t0, t [[,rtol], atol], res [,jac] [,hd])
24             [y, rd] = dae("root", initial, t0, t, res, ng, surface)
25             [y, rd [,hd]] = dae("root", initial, t0, t [[,rtol], atol], res [,jac], ng, surface [,hd])
26             [y, rd] = dae("root2", initial, t0, t, res, ng, surface)
27             [y, rd [,hd]] = dae("root2", initial, t0, t [[,rtol], atol], res [,jac], ng, surface [, psol, pjac] [, hd])
28         </synopsis>
29     </refsynopsisdiv>
30     <refsection>
31         <title>Paramètres</title>
32         <variablelist>
33             <varlistentry>
34                 <term>initial</term>
35                 <listitem>
36                     <para>
37                         vecteur colonne. Peut valoir <literal>x0</literal> ou
38                         <literal>[x0;xdot0]</literal>, où <literal>x0</literal> est la
39                         valeur de l'état au temps initial <literal>t0</literal> et
40                         <literal>xdot0</literal> est la valeur (ou une estimation)
41                         de l'état dérivé au temps initial (voir ci-après).
42                     </para>
43                 </listitem>
44             </varlistentry>
45             <varlistentry>
46                 <term>t0</term>
47                 <listitem>
48                     <para>réel, le temps initial.</para>
49                 </listitem>
50             </varlistentry>
51             <varlistentry>
52                 <term>t</term>
53                 <listitem>
54                     <para>scalaire ou vecteur réel. Les instants auxquels la solution est voulue.
55                         La solution peut s'obtenir à chaque étape de dae en initialisant
56                         <literal>
57                             <link linkend="daeoptions">%DAEOPTIONS</link>(2)=1
58                         </literal>
59                         .
60                     </para>
61                 </listitem>
62             </varlistentry>
63             <varlistentry>
64                 <term>atol</term>
65                 <listitem>
66                     <para>
67                         un scalaire ou vecteur colonne réel de même taille que
68                         <literal>x0</literal>, la tolérance absolue de la solution.
69                         Si <literal>atol</literal> est un vecteur, les tolérances
70                         sont spécifiées pour chaque composante de l'état.
71                     </para>
72                 </listitem>
73             </varlistentry>
74             <varlistentry>
75                 <term>rtol</term>
76                 <listitem>
77                     <para>
78                         un scalaire ou vecteur colonne réel de même taille que
79                         <literal>x0</literal>, la tolérance relative de la solution.
80                         Si <literal>atol</literal> est un vecteur, les tolérances
81                         sont spécifiées pour chaque composante de l'état.
82                     </para>
83                 </listitem>
84             </varlistentry>
85             <varlistentry>
86                 <term>res</term>
87                 <listitem>
88                     <para>
89                         une fonction <link linkend="external" role="" version="">external</link>, calcule la valeur de
90                         <literal>g(t,y,ydot)</literal>. Peut être
91                     </para>
92                     <variablelist>
93                         <varlistentry>
94                             <term>une fonction Scilab</term>
95                             <listitem>
96                                 <para>Dans ce cas, sa séquence d'appel doit être
97                                     <literal>[r,ires]=res(t,x,xdot)</literal> et
98                                     <literal>res</literal> doit retourner le résidu
99                                     <literal>r=g(t,x,xdot)</literal> et un drapeau d'erreur
100                                     <literal>ires</literal>.
101                                 </para>
102                                 <para>
103                                     <literal>ires = 0</literal> si
104                                     <literal>res</literal> est bien parvenu à calculer
105                                     <literal>r</literal>.
106                                 </para>
107                                 <para>
108                                     <literal>ires = -1</literal> si le résidu est
109                                     localement non-défini pour <literal>g(t,x,xdot)</literal>.
110                                 </para>
111                                 <para>
112                                     <literal>ires =-2</literal> si le paramètres ne sont pas admissibles.
113                                 </para>
114                             </listitem>
115                         </varlistentry>
116                         <varlistentry>
117                             <term>une liste</term>
118                             <listitem>
119                                 <para>
120                                     Cette forme d'<link linkend="external">external</link> sert à passer des paramètres à la fonction.
121                                     Elle doit se présenter comme suit:
122                                 </para>
123                                 <programlisting role="no-scilab-exec"><![CDATA[
124 list(res,p1,p2,...)
125  ]]></programlisting>
126                                 <para>où la séquence d'appel de la fonction
127                                     <literal>res</literal> est
128                                 </para>
129                                 <programlisting role="no-scilab-exec"><![CDATA[
130 r=res(t,y,ydot,p1,p2,...)
131  ]]></programlisting>
132                                 <para>
133                                     <literal>res</literal> retourne toujours le résidu comme fonction de
134                                     <literal>(t,x,xdot,x1,x2,...)</literal>, et
135                                     <literal>p1, p2,...</literal> sont des paramètres de la fonction.
136                                 </para>
137                             </listitem>
138                         </varlistentry>
139                         <varlistentry>
140                             <term>une chaîne de caractères</term>
141                             <listitem>
142                                 <para>
143                                     Elle doit se référer au nom d'une routine C ou Fortran.
144                                     Supposant que &lt;<literal>r_name</literal>&gt; est le nom donné,
145                                 </para>
146                                 <itemizedlist>
147                                     <listitem>
148                                         <para>
149                                             La séquence d'appel en Fortran doit être
150                                         </para>
151                                         <para>
152                                             <literal>&lt;r_name&gt;(t,x,xdot,res,ires,rpar,ipar)</literal>
153                                         </para>
154                                         <para>
155                                             <literal>double precision
156                                                 t,x(*),xdot(*),res(*),rpar(*)
157                                             </literal>
158                                         </para>
159                                         <para>
160                                             <literal>integer ires,ipar(*)</literal>
161                                         </para>
162                                     </listitem>
163                                     <listitem>
164                                         <para>
165                                             La séquence d'appel en C doit être
166                                         </para>
167                                         <para>
168                                             <literal>C2F(&lt;r_name&gt;)(double *t, double *x, double
169                                                 *xdot, double *res, integer *ires, double *rpar, integer
170                                                 *ipar)
171                                             </literal>
172                                         </para>
173                                     </listitem>
174                                 </itemizedlist>
175                                 <para>où</para>
176                                 <itemizedlist>
177                                     <listitem>
178                                         <para>
179                                             <literal>t</literal> est la valeur actuelle du temps
180                                         </para>
181                                     </listitem>
182                                     <listitem>
183                                         <para>
184                                             <literal>x</literal> est la valeur de l'état
185                                         </para>
186                                     </listitem>
187                                     <listitem>
188                                         <para>
189                                             <literal>xdot</literal> est la valeur de l'état dérivé
190                                         </para>
191                                     </listitem>
192                                     <listitem>
193                                         <para>
194                                             <literal>res</literal> la valeur du résidu
195                                         </para>
196                                     </listitem>
197                                     <listitem>
198                                         <para>
199                                             <literal>ires</literal> indicateur de complétion
200                                         </para>
201                                     </listitem>
202                                     <listitem>
203                                         <para>
204                                             <literal>rpar</literal> est un vecteur de paramètres nécéssaires
205                                             mais non initialisables par <literal>dae</literal>.
206                                         </para>
207                                     </listitem>
208                                     <listitem>
209                                         <para>
210                                             <literal>ipar</literal> est un vecteur de paramètres entiers
211                                             nécéssaires mais non initialisables par <literal>dae</literal>
212                                         </para>
213                                     </listitem>
214                                 </itemizedlist>
215                             </listitem>
216                         </varlistentry>
217                     </variablelist>
218                 </listitem>
219             </varlistentry>
220             <varlistentry>
221                 <term>jac</term>
222                 <listitem>
223                     <para>
224                         une fonction <link linkend="external">external</link>, calcule la valeur de
225                         <literal>dg/dx+cj*dg/dxdot</literal> pour un paramètre
226                         <literal>cj</literal> donné. Elle peut être:
227                     </para>
228                     <variablelist>
229                         <varlistentry>
230                             <term>une fonction Scilab</term>
231                             <listitem>
232                                 <para>Sa séquence d'appel doit être
233                                     <literal>r=jac(t,x,xdot,cj)</literal> et la fonction
234                                     <literal>jac</literal> doit retourner
235                                     <literal>r=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot</literal>
236                                     où <literal>cj</literal> est un scalaire réel.
237                                 </para>
238                             </listitem>
239                         </varlistentry>
240                         <varlistentry>
241                             <term>une liste</term>
242                             <listitem>
243                                 <para>
244                                     Cette forme d'<link linkend="external">external</link> est utilisée pour passer des paramètres à la fonction.
245                                     Elle doit se présenter comme suit :
246                                 </para>
247                                 <programlisting role="no-scilab-exec"><![CDATA[
248 list(jac,p1,p2,...)
249  ]]></programlisting>
250                                 <para>où la séquence d'appel de la fonction
251                                     <literal>jac</literal> est
252                                 </para>
253                                 <programlisting role="no-scilab-exec"><![CDATA[
254 r=jac(t,x,xdot,p1,p2,...)
255  ]]></programlisting>
256                                 <para>
257                                     <literal>jac</literal> retourne toujours
258                                     <literal>dg/dx+cj*dg/dxdot</literal> comme fonction de
259                                     <literal>(t,x,xdot,cj,p1,p2,...)</literal>.
260                                 </para>
261                             </listitem>
262                         </varlistentry>
263                         <varlistentry>
264                             <term>une chaîne de caractères</term>
265                             <listitem>
266                                 <para>elle doit se référer au nom d'une routine C ou Fortran.
267                                     Supposant que &lt;j_name&gt; est le nom donné,
268                                 </para>
269                                 <itemizedlist>
270                                     <listitem>
271                                         <para>
272                                             La séquence d'appel en Fortran doit être
273                                         </para>
274                                         <para>
275                                             <literal>&lt;j_name&gt;(t, x, xdot, r, cj, ires,
276                                                 rpar, ipar)
277                                             </literal>
278                                         </para>
279                                         <para>
280                                             double precision <literal>t, x(*), xdot(*), r(*),
281                                                 ci, rpar(*)
282                                             </literal>
283                                         </para>
284                                         <para>
285                                             integer <literal>ires, ipar(*)</literal>
286                                         </para>
287                                     </listitem>
288                                     <listitem>
289                                         <para>
290                                             La séquence d'appel en C doit être
291                                         </para>
292                                         <para>
293                                             <literal>C2F(&lt;j_name&gt;)(double *t, double *x, double
294                                                 *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)
295                                             </literal>
296                                         </para>
297                                     </listitem>
298                                 </itemizedlist>
299                                 <para>
300                                     où <literal>t, x, xdot, ires, rpar, ipar</literal>
301                                     ont la même définition qu'au-dessus, <literal>r</literal>
302                                     est la valeur du résultat
303                                 </para>
304                             </listitem>
305                         </varlistentry>
306                     </variablelist>
307                 </listitem>
308             </varlistentry>
309             <varlistentry>
310                 <term>surface</term>
311                 <listitem>
312                     <para>
313                         une fonction <link linkend="external">external</link>, calcule a valeur du vecteur colonne
314                         <literal>surface(t,x)</literal> à <literal>ng</literal> éléments.
315                         Chaque élément définit une surface. Elle peut être
316                     </para>
317                     <variablelist>
318                         <varlistentry>
319                             <term>une fonction Scilab</term>
320                             <listitem>
321                                 <para>Sa séquence d'appel doit être
322                                     <literal>r=surface(t,x)</literal>, cette fonction doit retourner
323                                     un vecteur colonne à <literal>ng</literal> éléments.
324                                 </para>
325                             </listitem>
326                         </varlistentry>
327                         <varlistentry>
328                             <term>une liste</term>
329                             <listitem>
330                                 <para>
331                                     Cette forme d'<link linkend="external">external</link> est utilisé pour
332                                     passer des paramètres à la fonction.
333                                     Elle doit se présenter comme suit :
334                                 </para>
335                                 <programlisting role="no-scilab-exec"><![CDATA[
336 list(surface,p1,p2,...)
337  ]]></programlisting>
338                                 <para>où la qéquence d'appel de la fonction
339                                     <literal>surface</literal> est
340                                 </para>
341                                 <programlisting role="no-scilab-exec"><![CDATA[
342 r=surface(t,x,p1,p2,...)
343  ]]></programlisting>
344                             </listitem>
345                         </varlistentry>
346                         <varlistentry>
347                             <term>une chaîne de caractères</term>
348                             <listitem>
349                                 <para>elle doit se référer au nom d'une routine C ou Fortran.
350                                     Supposant que &lt;s_name&gt; est le nom donné,
351                                 </para>
352                                 <itemizedlist>
353                                     <listitem>
354                                         <para>
355                                             La séquence d'appel en Fortran doit être
356                                         </para>
357                                         <para>
358                                             <literal>&lt;s_name&gt;(nx, t, x, ng, r, rpar,
359                                                 ipar)
360                                             </literal>
361                                         </para>
362                                         <para>
363                                             <literal>double precision t, x(*), r(*),
364                                                 rpar(*)
365                                             </literal>
366                                         </para>
367                                         <para>
368                                             <literal>integer nx, ng,ipar(*)</literal>
369                                         </para>
370                                     </listitem>
371                                     <listitem>
372                                         <para>
373                                             La séquence d'appel en C doit être
374                                         </para>
375                                         <para>
376                                             <literal>C2F(&lt;s_name&gt;)(double *t, double *x, double
377                                                 *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)
378                                             </literal>
379                                         </para>
380                                     </listitem>
381                                 </itemizedlist>
382                                 <para>
383                                     où <literal>t, x, rpar, ipar</literal> ont les mêmes définitions qu'au-dessus,
384                                     <literal>ng</literal> est le nombre de surfaces,
385                                     <literal>nx</literal> est la dimension de l'état et
386                                     <literal>r</literal> est le vecteur résultat.
387                                 </para>
388                             </listitem>
389                         </varlistentry>
390                     </variablelist>
391                 </listitem>
392             </varlistentry>
393             <varlistentry>
394                 <term>rd</term>
395                 <listitem>
396                     <para>
397                         un vecteur à deux entrées <literal>[times num]</literal> où
398                         <literal>times</literal> est la valeur des temps auquel les surfaces sont traversées,
399                         <literal>num</literal> est le nombre de surfaces traversées.
400                     </para>
401                 </listitem>
402             </varlistentry>
403             <varlistentry>
404                 <term>pjac</term>
405                 <listitem>
406                     <para>
407                         external (fonction, liste or chaîne de caractères). Voir <link linkend="daskr">daskr</link>
408                     </para>
409                 </listitem>
410             </varlistentry>
411             <varlistentry>
412                 <term>psol</term>
413                 <listitem>
414                     <para>
415                         external (fonction, liste or chaîne de caractères). Voir <link linkend="daskr">daskr</link>
416                     </para>
417                 </listitem>
418             </varlistentry>
419             <varlistentry>
420                 <term>hd</term>
421                 <listitem>
422                     <para>vecteur réel, stocke en sortie le contexte de
423                         <literal>dae</literal>. Peut être utilisé comme paramètre d'entrée
424                         pour reprendre l'intégration (reprise à chaud).
425                     </para>
426                 </listitem>
427             </varlistentry>
428             <varlistentry>
429                 <term>y</term>
430                 <listitem>
431                     <para>
432                         matrice réelle. Si
433                         <literal>
434                             <link linkend="daeoptions">%DAEOPTIONS</link>(2)=1
435                         </literal>
436                         ,
437                         chaque colonne est le vecteur <literal>[t;x(t);xdot(t)]</literal> où
438                         <literal>t</literal> est le temps auquel la solution à été calculée.
439                         Sinon <literal>y</literal> est le vecteur
440                         <literal>[x(t);xdot(t)]</literal>.
441                     </para>
442                 </listitem>
443             </varlistentry>
444         </variablelist>
445     </refsection>
446     <refsection>
447         <title>Description</title>
448         <para>
449             La fonction <literal>dae</literal> est une gateway construite sur les solveurs
450             <link linkend="dassl">dassl</link>, <link linkend="dasrt">dasrt</link>
451             et <link linkend="daskr">daskr</link>
452             destinés à l'intégration implicite d'équations différentielles.
453         </para>
454         <para>
455             L'option <literal>"root"</literal> appelle la routine <link linkend="dasrt">dasrt</link>,
456             et <literal>"root2"</literal> appelle <link linkend="dasrt">daskr</link>.
457         </para>
458         <programlisting role="no-scilab-exec"><![CDATA[
459 g(t, x, xdot) = 0
460 x(t0) = x0 and xdot(t0) = xdot0
461  ]]></programlisting>
462         <para>
463             Si <literal>xdot0</literal> n'est pas donné en paramètre initial,
464             <literal>dae</literal> tente de la calculer en résolvant
465             <literal>g(t,x0,xdot0)=0</literal>.
466         </para>
467         <para>
468             Si <literal>xdot0</literal> est donné en paramètre initial, il peut soit
469             satifaire <literal>g(t,x0,xdot0)=0</literal>, soit en être une approximation.
470             Dans le dernier cas,
471             <link linkend="daeoptions">%DAEOPTIONS</link>(7) doit être initialisé à 1.
472         </para>
473         <para>Des exemples détaillés utilisant des externals codées en C et en Scilab se trouvent dans
474             <literal>modules/differential_equations/tests/unit_tests/dassldasrt.tst</literal> et
475             <literal>modules/differential_equations/tests/unit_tests/daskr.tst</literal>
476         </para>
477     </refsection>
478     <refsection>
479         <title>Exemples</title>
480         <para>
481             Exemple #1: dassl (pas de traversée de surface)
482         </para>
483         <programlisting role="example"><![CDATA[
484 // Exemple avec du code Scilab
485 --------------------------------------------------
486 function [r, ires] = chemres(t, y, yd)
487     r(1) = -0.04*y(1) + 1d4*y(2)*y(3) - yd(1);
488     r(2) =  0.04*y(1) - 1d4*y(2)*y(3) - 3d7*y(2)*y(2) - yd(2);
489     r(3) =       y(1) +     y(2)      + y(3)-1;
490     ires =  0;
491 endfunction
492
493 function pd = chemjac(x, y, yd, cj)
494     pd = [-0.04-cj , 1d4*y(3)               , 1d4*y(2);
495            0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
496            1       , 1                      , 1       ]
497 endfunction
498
499 x0 = [1; 0; 0];
500 xd0 = [-0.04; 0.04; 0];
501 t = [1.d-5:0.02:.4, 0.41:.1:4, 40, 400, 4000, 40000, 4d5, 4d6, 4d7, 4d8, 4d9, 4d10];
502
503 y = dae([x0, xd0], 0, t, chemres); // Retourne les points d'observation demandés
504
505 %DAEOPTIONS = list([], 1, [], [], [], 0, 0); // Demande à dae les points à retourner
506 y = dae([x0, xd0], 0, 4d10, chemres); // Sans jacobian
507 y = dae([x0, xd0], 0, 4d10, chemres, chemjac); // Avec jacobien
508  ]]></programlisting>
509         <para>
510             Exemple #2: dasrt ("root")
511         </para>
512         <programlisting role="example"><![CDATA[
513 // Exemple avec du code C (compilateur C requis)
514 --------------------------------------------------
515 bOK = haveacompiler();
516 if bOK <> %t
517     [btn] = messagebox(["Vous avez besoin nd'un compilateur C pour cet exemple."; "Arrêt de l'exécution."], "Problème de Software", 'info');
518     return
519 end
520
521 //-1- Crée les codes C dans TMPDIR - équation de Vanderpol, forme implicite
522 code = ['#include <math.h>'
523       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
524       '{res[0] = yd[0] - y[1];'
525       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
526       ' '
527       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
528       '{pd[0] = *cj - 0.0;'
529       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
530       ' pd[2] =     - 1.0;'
531       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
532       ' '
533       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
534       '{ groot[0] = y[0];}']
535 previous_dir = pwd();
536 cd TMPDIR;
537 mputl(code, 't22.c')
538
539 //-2- Compile et charge
540 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
541 exec('t22loader.sce')
542
543 //-3- Exécute
544 rtol = [1.d-6; 1.d-6];
545 atol = [1.d-6; 1.d-4];
546 t0 = 0; t = [20:20:200];
547 y0 = [2; 0]; y0d = [0; -2];
548 ng = 1;
549
550 // Simulation simple
551 t = 0:0.003:300;
552 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
553 clf(); plot(yy(1, :), yy(2, :))
554 // Trouve le premier point où yy(1) = 0
555 [yy, nn, hotd] = dae("root", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22');
556 plot(yy(1, 1), yy(2, 1), 'r+')
557 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
558
559 // Reprise à chaud pour le point suivant
560 t01 = nn(1);
561 [pp, qq] = size(yy);
562 y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
563 [yy, nn, hotd] = dae("root", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', hotd);
564 plot(yy(1, 1), yy(2, 1), 'r+')
565 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
566 cd(previous_dir);
567  ]]></programlisting>
568         <scilab:image><![CDATA[
569 code = ['#include <math.h>'
570       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
571       '{res[0] = yd[0] - y[1];'
572       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
573       ' '
574       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
575       '{pd[0] = *cj - 0.0;'
576       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
577       ' pd[2] =     - 1.0;'
578       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
579       ' '
580       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
581       '{ groot[0] = y[0];}']
582 previous_dir = pwd();
583 cd TMPDIR;
584 mputl(code, 't22.c')
585 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
586 exec('t22loader.sce')
587 rtol = [1.d-6; 1.d-6];
588 atol = [1.d-6; 1.d-4];
589 t0 = 0; t = [20:20:200];
590 y0 = [2; 0]; y0d = [0; -2];
591 ng = 1;
592 t = 0:0.003:300;
593 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
594 clf(); plot(yy(1, :), yy(2, :))
595 [yy, nn, hotd] = dae("root", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22');
596 plot(yy(1, 1), yy(2, 1), 'r+')
597 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
598 t01 = nn(1);
599 [pp, qq] = size(yy);
600 y01 = yy(2:3, qq);
601 y0d1 = yy(3:4, qq);
602 [yy, nn, hotd] = dae("root", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', hotd);
603 plot(yy(1, 1), yy(2, 1), 'r+')
604 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
605 cd(previous_dir);
606  ]]></scilab:image>
607         <para>
608             Exemple #3: daskr ("root2"), avec les routines 'psol' et 'pjac' par défaut
609         </para>
610         <programlisting role="example"><![CDATA[
611 // Exemple avec du code C (compilateur C requis)
612 --------------------------------------------------
613 bOK = haveacompiler();
614 if bOK <> %t
615     [btn] = messagebox(["Vous avez besoin nd'un compilateur C pour cet exemple."; "Arrêt de l'exécution."], "Problème de Software", 'info');
616     return
617 end
618
619 //-1- Crée les codes C dans TMPDIR - équation de Vanderpol, forme implicite
620 code = ['#include <math.h>'
621       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
622       '{res[0] = yd[0] - y[1];'
623       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
624       ' '
625       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
626       '{pd[0] = *cj - 0.0;'
627       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
628       ' pd[2] =     - 1.0;'
629       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
630       ' '
631       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
632       '{ groot[0] = y[0];}']
633 previous_dir = pwd();
634 cd TMPDIR;
635 mputl(code, 't22.c')
636
637 //-2- Compile et charge
638 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
639 exec('t22loader.sce')
640
641 //-3- Exécute
642 rtol = [1.d-6; 1.d-6];
643 atol = [1.d-6; 1.d-4];
644 t0 = 0; t = [20:20:200];
645 y0 = [2; 0]; y0d = [0; -2];
646 ng = 1;
647
648 // Simulation simple
649 t = 0:0.003:300;
650 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
651 clf(); plot(yy(1, :), yy(2, :))
652 // Trouve le premier point où yy(1) = 0
653 %DAEOPTIONS = list([] , 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
654 [yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
655 plot(yy(1, 1), yy(2, 1), 'r+')
656 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
657
658 // Reprise à chaud pour le point suivant
659 t01 = nn(1);
660 [pp, qq] = size(yy);
661 y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
662 [yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
663 plot(yy(1, 1), yy(2, 1), 'r+')
664 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
665 cd(previous_dir);
666  ]]></programlisting>
667         <scilab:image><![CDATA[
668 code = ['#include <math.h>'
669       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
670       '{res[0] = yd[0] - y[1];'
671       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
672       ' '
673       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
674       '{pd[0] = *cj - 0.0;'
675       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
676       ' pd[2] =     - 1.0;'
677       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
678       ' '
679       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
680       '{ groot[0] = y[0];}']
681 previous_dir = pwd();
682 cd TMPDIR;
683 mputl(code, 't22.c')
684 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
685 exec('t22loader.sce')
686 rtol = [1.d-6; 1.d-6];
687 atol = [1.d-6; 1.d-4];
688 t0 = 0; t = [20:20:200];
689 y0 = [2; 0]; y0d = [0; -2];
690 ng = 1;
691 t = 0:0.003:300;
692 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
693 clf(); plot(yy(1, :), yy(2, :))
694 %DAEOPTIONS = list([], 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
695 [yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
696 plot(yy(1, 1), yy(2, 1), 'r+')
697 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
698 t01 = nn(1);
699 [pp, qq] = size(yy);
700 y01 = yy(2:3, qq);
701 y0d1 = yy(3:4, qq);
702 [yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
703 plot(yy(1, 1), yy(2, 1), 'r+')
704 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
705 cd(previous_dir);
706  ]]></scilab:image>
707     </refsection>
708     <refsection role="see also">
709         <title>Voir aussi</title>
710         <simplelist type="inline">
711             <member>
712                 <link linkend="ode">ode</link>
713             </member>
714             <member>
715                 <link linkend="daeoptions">daeoptions</link>
716             </member>
717             <member>
718                 <link linkend="dassl">dassl</link>
719             </member>
720             <member>
721                 <link linkend="dasrt">dasrt</link>
722             </member>
723             <member>
724                 <link linkend="daskr">daskr</link>
725             </member>
726             <member>
727                 <link linkend="impl">impl</link>
728             </member>
729             <member>
730                 <link linkend="fort">fort</link>
731             </member>
732             <member>
733                 <link linkend="link">link</link>
734             </member>
735             <member>
736                 <link linkend="external">external</link>
737             </member>
738         </simplelist>
739     </refsection>
740 </refentry>