Differential_equations help: fix dae tolerances
[scilab.git] / scilab / modules / differential_equations / help / en_US / 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>Differential algebraic equations solver</refpurpose>
19     </refnamediv>
20     <refsynopsisdiv>
21         <title>Calling Sequence</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>Arguments</title>
32         <variablelist>
33             <varlistentry>
34                 <term>initial</term>
35                 <listitem>
36                     <para>
37                         a column vector. It may be equal to <literal>x0</literal> or
38                         <literal>[x0;xdot0]</literal>. Where <literal>x0</literal> is the
39                         state value at initial time <literal>t0</literal> and
40                         <literal>xdot0</literal> is the initial state derivative value or an
41                         estimation of it (see below).
42                     </para>
43                 </listitem>
44             </varlistentry>
45             <varlistentry>
46                 <term>t0</term>
47                 <listitem>
48                     <para>a real number, the initial time.</para>
49                 </listitem>
50             </varlistentry>
51             <varlistentry>
52                 <term>t</term>
53                 <listitem>
54                     <para>a real scalar or vector. Gives instants for which you want the
55                         solution. Note that you can get solution at each dae's step point by
56                         setting
57                         <literal>
58                             <link linkend="daeoptions">%DAEOPTIONS</link>(2)=1
59                         </literal>
60                         .
61                     </para>
62                 </listitem>
63             </varlistentry>
64             <varlistentry>
65                 <term>atol</term>
66                 <listitem>
67                     <para>a real scalar or a column vector of same size as
68                         <literal>x0</literal>, the absolute error
69                         tolerance of solution. If <literal>atol</literal> is a vector the
70                         tolerances are specified for each component of the state.
71                     </para>
72                 </listitem>
73             </varlistentry>
74             <varlistentry>
75                 <term>rtol</term>
76                 <listitem>
77                     <para>
78                         a real scalar or a column vector of same size as
79                         <literal>x0</literal>, the relative error
80                         tolerance of solution. If <literal>rtol</literal> is a vector the
81                         tolerances are specified for each component of the state.
82                     </para>
83                 </listitem>
84             </varlistentry>
85             <varlistentry>
86                 <term>res</term>
87                 <listitem>
88                     <para>
89                         an <link linkend="external" role="" version="">external</link> function, computes the value of
90                         <literal>g(t,y,ydot)</literal>. It may be
91                     </para>
92                     <variablelist>
93                         <varlistentry>
94                             <term>a Scilab function</term>
95                             <listitem>
96                                 <para>In this case, its calling sequence must be
97                                     <literal>[r,ires]=res(t,x,xdot)</literal> and
98                                     <literal>res</literal> must return the residue
99                                     <literal>r=g(t,x,xdot)</literal> and error flag
100                                     <literal>ires</literal>.
101                                 </para>
102                                 <para>
103                                     <literal>ires = 0</literal> if
104                                     <literal>res</literal> succeeds to compute
105                                     <literal>r</literal>.
106                                 </para>
107                                 <para>
108                                     <literal>ires = -1</literal> if residue is
109                                     locally not defined for <literal>g(t,x,xdot)</literal>.
110                                 </para>
111                                 <para>
112                                     <literal>ires =-2</literal> if parameters are out of admissible
113                                     range.
114                                 </para>
115                             </listitem>
116                         </varlistentry>
117                         <varlistentry>
118                             <term>a list</term>
119                             <listitem>
120                                 <para>This form of external is used to pass parameters to the
121                                     function. It must be as follows:
122                                 </para>
123                                 <programlisting role="no-scilab-exec"><![CDATA[
124 list(res,p1,p2,...)
125  ]]></programlisting>
126                                 <para>where the calling sequence of the function
127                                     <literal>res</literal> is now
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> still returns the residual value
134                                     as a function of <literal>(t,x,xdot,x1,x2,...)</literal>, and
135                                     <literal>p1, p2,...</literal> are function parameters.
136                                 </para>
137                             </listitem>
138                         </varlistentry>
139                         <varlistentry>
140                             <term>a character string</term>
141                             <listitem>
142                                 <para>
143                                     it must refer to the name of a C or Fortran routine.
144                                     Assuming that &lt;<literal>r_name</literal>&gt; is the given name,
145                                 </para>
146                                 <itemizedlist>
147                                     <listitem>
148                                         <para>
149                                             The Fortran calling sequence must be
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                                             The C calling sequence must be
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>where</para>
176                                 <itemizedlist>
177                                     <listitem>
178                                         <para>
179                                             <literal>t</literal> is the current time
180                                             value
181                                         </para>
182                                     </listitem>
183                                     <listitem>
184                                         <para>
185                                             <literal>x</literal> the state array
186                                         </para>
187                                     </listitem>
188                                     <listitem>
189                                         <para>
190                                             <literal>xdot</literal> the array of state
191                                             derivatives
192                                         </para>
193                                     </listitem>
194                                     <listitem>
195                                         <para>
196                                             <literal>res</literal> the array of residuals
197                                         </para>
198                                     </listitem>
199                                     <listitem>
200                                         <para>
201                                             <literal>ires</literal> the execution
202                                             indicator
203                                         </para>
204                                     </listitem>
205                                     <listitem>
206                                         <para>
207                                             <literal>rpar</literal> is the array of floating
208                                             point parameter values, needed but cannot be set by the
209                                             <literal>dae</literal> function
210                                         </para>
211                                     </listitem>
212                                     <listitem>
213                                         <para>
214                                             <literal>ipar</literal> is the array of floating
215                                             integer parameter values, needed but cannot be set by the
216                                             <literal>dae</literal> function
217                                         </para>
218                                     </listitem>
219                                 </itemizedlist>
220                             </listitem>
221                         </varlistentry>
222                     </variablelist>
223                 </listitem>
224             </varlistentry>
225             <varlistentry>
226                 <term>jac</term>
227                 <listitem>
228                     <para>
229                         an <link linkend="external">external</link>, computes the value of
230                         <literal>dg/dx+cj*dg/dxdot</literal> for a given value of parameter
231                         <literal>cj</literal>. It may be
232                     </para>
233                     <variablelist>
234                         <varlistentry>
235                             <term>a Scilab function</term>
236                             <listitem>
237                                 <para>Its calling sequence must be
238                                     <literal>r=jac(t,x,xdot,cj)</literal> and the
239                                     <literal>jac</literal> function must return
240                                     <literal>r=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot</literal>
241                                     where <literal>cj</literal> is a real scalar.
242                                 </para>
243                             </listitem>
244                         </varlistentry>
245                         <varlistentry>
246                             <term>a list</term>
247                             <listitem>
248                                 <para>This form of external is used to pass parameters to the
249                                     function. It must be as follows:
250                                 </para>
251                                 <programlisting role="no-scilab-exec"><![CDATA[
252 list(jac,p1,p2,...)
253  ]]></programlisting>
254                                 <para>where the calling sequence of the function
255                                     <literal>jac</literal> is now
256                                 </para>
257                                 <programlisting role="no-scilab-exec"><![CDATA[
258 r=jac(t,x,xdot,p1,p2,...)
259  ]]></programlisting>
260                                 <para>
261                                     <literal>jac</literal> still returns
262                                     <literal>dg/dx+cj*dg/dxdot</literal> as a function of
263                                     <literal>(t,x,xdot,cj,p1,p2,...)</literal>.
264                                 </para>
265                             </listitem>
266                         </varlistentry>
267                         <varlistentry>
268                             <term>a character string</term>
269                             <listitem>
270                                 <para>it must refer to the name of a C or Fortran routine.
271                                     Assuming that &lt;j_name&gt; is the given name,
272                                 </para>
273                                 <itemizedlist>
274                                     <listitem>
275                                         <para>
276                                             The Fortran calling sequence must be
277                                         </para>
278                                         <para>
279                                             <literal>&lt;j_name&gt;(t, x, xdot, r, cj, ires,
280                                                 rpar, ipar)
281                                             </literal>
282                                         </para>
283                                         <para>
284                                             double precision <literal>t, x(*), xdot(*), r(*),
285                                                 ci, rpar(*)
286                                             </literal>
287                                         </para>
288                                         <para>
289                                             integer <literal>ires, ipar(*)</literal>
290                                         </para>
291                                     </listitem>
292                                     <listitem>
293                                         <para>
294                                             The C calling sequence must be
295                                         </para>
296                                         <para>
297                                             <literal>C2F(&lt;j_name&gt;)(double *t, double *x, double
298                                                 *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)
299                                             </literal>
300                                         </para>
301                                     </listitem>
302                                 </itemizedlist>
303                                 <para>
304                                     where <literal>t, x, xdot, ires, rpar, ipar</literal>
305                                     have similar definition as above, <literal>r</literal> is the results
306                                     array
307                                 </para>
308                             </listitem>
309                         </varlistentry>
310                     </variablelist>
311                 </listitem>
312             </varlistentry>
313             <varlistentry>
314                 <term>surface</term>
315                 <listitem>
316                     <para>
317                         an <link linkend="external">external</link>, computes the value of the column vector
318                         <literal>surface(t,x)</literal> with <literal>ng</literal>
319                         components. Each component defines a surface.
320                     </para>
321                     <variablelist>
322                         <varlistentry>
323                             <term>a Scilab function</term>
324                             <listitem>
325                                 <para>Its calling sequence must be
326                                     <literal>r=surface(t,x)</literal>, this function must return a
327                                     vector with <literal>ng</literal> elements.
328                                 </para>
329                             </listitem>
330                         </varlistentry>
331                         <varlistentry>
332                             <term>a list</term>
333                             <listitem>
334                                 <para>
335                                     This form of <link linkend="external">external</link> is used to pass parameters to the
336                                     function. It must be as follows:
337                                 </para>
338                                 <programlisting role="no-scilab-exec"><![CDATA[
339 list(surface,p1,p2,...)
340  ]]></programlisting>
341                                 <para>where the calling sequence of the function
342                                     <literal>surface</literal> is now
343                                 </para>
344                                 <programlisting role="no-scilab-exec"><![CDATA[
345 r=surface(t,x,p1,p2,...)
346  ]]></programlisting>
347                             </listitem>
348                         </varlistentry>
349                         <varlistentry>
350                             <term>a character string</term>
351                             <listitem>
352                                 <para>it must refer to the name of a C or Fortran routine.
353                                     Assuming that &lt;s_name&gt; is the given name,
354                                 </para>
355                                 <itemizedlist>
356                                     <listitem>
357                                         <para>
358                                             The Fortran calling sequence must be
359                                         </para>
360                                         <para>
361                                             <literal>&lt;s_name&gt;(nx, t, x, ng, r, rpar,
362                                                 ipar)
363                                             </literal>
364                                         </para>
365                                         <para>
366                                             <literal>double precision t, x(*), r(*),
367                                                 rpar(*)
368                                             </literal>
369                                         </para>
370                                         <para>
371                                             <literal>integer nx, ng,ipar(*)</literal>
372                                         </para>
373                                     </listitem>
374                                     <listitem>
375                                         <para>
376                                             The C calling sequence must be
377                                         </para>
378                                         <para>
379                                             <literal>C2F(&lt;s_name&gt;)(double *t, double *x, double
380                                                 *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)
381                                             </literal>
382                                         </para>
383                                     </listitem>
384                                 </itemizedlist>
385                                 <para>
386                                     where <literal>t, x, rpar, ipar</literal> have similar
387                                     definition as above, <literal>ng</literal> is the number of
388                                     surfaces, <literal>nx</literal> the dimension of the state and
389                                     <literal>r</literal> is the results array.
390                                 </para>
391                             </listitem>
392                         </varlistentry>
393                     </variablelist>
394                 </listitem>
395             </varlistentry>
396             <varlistentry>
397                 <term>rd</term>
398                 <listitem>
399                     <para>
400                         a vector with two entries <literal>[times num]</literal> where
401                         <literal>times</literal> is the value of the time at which the
402                         surface is crossed, <literal>num</literal> is the number of the
403                         crossed surface
404                     </para>
405                 </listitem>
406             </varlistentry>
407             <varlistentry>
408                 <term>pjac</term>
409                 <listitem>
410                     <para>
411                         external (function, list or string). See <link linkend="daskr">daskr</link>
412                     </para>
413                 </listitem>
414             </varlistentry>
415             <varlistentry>
416                 <term>psol</term>
417                 <listitem>
418                     <para>
419                         external (function, list or string). See <link linkend="daskr">daskr</link>
420                     </para>
421                 </listitem>
422             </varlistentry>
423             <varlistentry>
424                 <term>hd</term>
425                 <listitem>
426                     <para>a real vector, as an output it stores the
427                         <literal>dae</literal> context. It can be used as an input argument
428                         to resume integration (hot restart).
429                     </para>
430                 </listitem>
431             </varlistentry>
432             <varlistentry>
433                 <term>y</term>
434                 <listitem>
435                     <para>
436                         a real matrix. If
437                         <literal>
438                             <link linkend="daeoptions">%DAEOPTIONS</link>(2)=1
439                         </literal>
440                         ,
441                         each column is the vector <literal>[t;x(t);xdot(t)]</literal> where
442                         <literal>t</literal> is time index for which the solution has been
443                         computed. Else <literal>y</literal> is the vector
444                         <literal>[x(t);xdot(t)]</literal>.
445                     </para>
446                 </listitem>
447             </varlistentry>
448         </variablelist>
449     </refsection>
450     <refsection>
451         <title>Description</title>
452         <para>
453             The <literal>dae</literal> function is a gateway built above the
454             <link linkend="dassl">dassl</link>, <link linkend="dasrt">dasrt</link>
455             and <link linkend="daskr">daskr</link>
456             functions designed for implicit differential equations integration.
457         </para>
458         <para>
459             Option <literal>"root"</literal> calls the <link linkend="dasrt">dasrt</link> routine,
460             and <literal>"root2"</literal> calls <link linkend="dasrt">daskr</link>.
461         </para>
462         <programlisting role="no-scilab-exec"><![CDATA[
463 g(t, x, xdot) = 0
464 x(t0) = x0 and xdot(t0) = xdot0
465  ]]></programlisting>
466         <para>
467             If <literal>xdot0</literal> is not given in the <emphasis>initial</emphasis>
468             argument, the <literal>dae</literal> function tries to compute it solving
469             <literal>g(t,x0,xdot0)=0</literal>.
470         </para>
471         <para>
472             If <literal>xdot0</literal> is given in the <emphasis>initial</emphasis>
473             argument it may be either a compatible derivative
474             satisfying <literal>g(t,x0,xdot0)=0</literal> or an approximate value. In the latter case
475             <link linkend="daeoptions">%DAEOPTIONS</link>(7) must be set to 1.
476         </para>
477         <para>Detailed examples using Scilab and C coded externals are given in
478             <literal>modules/differential_equations/tests/unit_tests/dassldasrt.tst</literal> and
479             <literal>modules/differential_equations/tests/unit_tests/daskr.tst</literal>
480         </para>
481     </refsection>
482     <refsection>
483         <title>Examples</title>
484         <para>
485             Example #1: dassl (no root finding)
486         </para>
487         <programlisting role="example"><![CDATA[
488 // Example with Scilab code
489 --------------------------------------------------
490 function [r, ires] = chemres(t, y, yd)
491     r(1) = -0.04*y(1) + 1d4*y(2)*y(3) - yd(1);
492     r(2) =  0.04*y(1) - 1d4*y(2)*y(3) - 3d7*y(2)*y(2) - yd(2);
493     r(3) =       y(1) +     y(2)      + y(3)-1;
494     ires =  0;
495 endfunction
496
497 function pd = chemjac(x, y, yd, cj)
498     pd = [-0.04-cj , 1d4*y(3)               , 1d4*y(2);
499            0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
500            1       , 1                      , 1       ]
501 endfunction
502
503 x0 = [1; 0; 0];
504 xd0 = [-0.04; 0.04; 0];
505 t = [1.d-5:0.02:.4, 0.41:.1:4, 40, 400, 4000, 40000, 4d5, 4d6, 4d7, 4d8, 4d9, 4d10];
506
507 y = dae([x0, xd0], 0, t, chemres); // Returns requested observation time points
508
509 %DAEOPTIONS = list([], 1, [], [], [], 0, 0); // Ask  dae mesh points to be returned
510 y = dae([x0, xd0], 0, 4d10, chemres); // Without jacobian
511 y = dae([x0, xd0], 0, 4d10, chemres, chemjac); // With jacobian
512  ]]></programlisting>
513         <para>
514             Example #2: dasrt ("root")
515         </para>
516         <programlisting role="example"><![CDATA[
517 // Example with C code (C compiler needed)
518 --------------------------------------------------
519 bOK = haveacompiler();
520 if bOK <> %t
521     [btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info');
522     return
523 end
524
525 //-1- Create the C codes in TMPDIR - Vanderpol equation, implicit form
526 code = ['#include <math.h>'
527       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
528       '{res[0] = yd[0] - y[1];'
529       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
530       ' '
531       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
532       '{pd[0] = *cj - 0.0;'
533       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
534       ' pd[2] =     - 1.0;'
535       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
536       ' '
537       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
538       '{ groot[0] = y[0];}']
539 previous_dir = pwd();
540 cd TMPDIR;
541 mputl(code, 't22.c')
542
543 //-2- Compile and load them
544 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
545 exec('t22loader.sce')
546
547 //-3- Run
548 rtol = [1.d-6; 1.d-6];
549 atol = [1.d-6; 1.d-4];
550 t0 = 0; t = [20:20:200];
551 y0 = [2; 0]; y0d = [0; -2];
552 ng = 1;
553
554 // Simple simulation
555 t = 0:0.003:300;
556 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
557 clf(); plot(yy(1, :), yy(2, :))
558 // Find first point where yy(1) = 0
559 [yy, nn, hotd] = dae("root", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22');
560 plot(yy(1, 1), yy(2, 1), 'r+')
561 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
562
563 // Hot restart for next point
564 t01 = nn(1);
565 [pp, qq] = size(yy);
566 y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
567 [yy, nn, hotd] = dae("root", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', hotd);
568 plot(yy(1, 1), yy(2, 1), 'r+')
569 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
570 cd(previous_dir);
571  ]]></programlisting>
572         <scilab:image><![CDATA[
573 code = ['#include <math.h>'
574       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
575       '{res[0] = yd[0] - y[1];'
576       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
577       ' '
578       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
579       '{pd[0] = *cj - 0.0;'
580       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
581       ' pd[2] =     - 1.0;'
582       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
583       ' '
584       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
585       '{ groot[0] = y[0];}']
586 previous_dir = pwd();
587 cd TMPDIR;
588 mputl(code, 't22.c')
589 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
590 exec('t22loader.sce')
591 rtol = [1.d-6; 1.d-6];
592 atol = [1.d-6; 1.d-4];
593 t0 = 0; t = [20:20:200];
594 y0 = [2; 0]; y0d = [0; -2];
595 ng = 1;
596 t = 0:0.003:300;
597 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
598 clf(); plot(yy(1, :), yy(2, :))
599 [yy, nn, hotd] = dae("root", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22');
600 plot(yy(1, 1), yy(2, 1), 'r+')
601 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
602 t01 = nn(1);
603 [pp, qq] = size(yy);
604 y01 = yy(2:3, qq);
605 y0d1 = yy(3:4, qq);
606 [yy, nn, hotd] = dae("root", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', hotd);
607 plot(yy(1, 1), yy(2, 1), 'r+')
608 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
609 cd(previous_dir);
610  ]]></scilab:image>
611         <para>
612             Example #3: daskr ("root2"), using default 'psol' and 'pjac' routines
613         </para>
614         <programlisting role="example"><![CDATA[
615 // Example with C code (C compiler needed)
616 --------------------------------------------------
617 bOK = haveacompiler();
618 if bOK <> %t
619     [btn] = messagebox(["You need a C compiler for this example."; "Execution of this example is canceled."], "Software problem", 'info');
620     return
621 end
622
623 //-1- Create the C codes in TMPDIR - Vanderpol equation, implicit form
624 code = ['#include <math.h>'
625       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
626       '{res[0] = yd[0] - y[1];'
627       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
628       ' '
629       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
630       '{pd[0] = *cj - 0.0;'
631       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
632       ' pd[2] =     - 1.0;'
633       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
634       ' '
635       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
636       '{ groot[0] = y[0];}']
637 previous_dir = pwd();
638 cd TMPDIR;
639 mputl(code, 't22.c')
640
641 //-2- Compile and load them
642 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
643 exec('t22loader.sce')
644
645 //-3- Run
646 rtol = [1.d-6; 1.d-6];
647 atol = [1.d-6; 1.d-4];
648 t0 = 0; t = [20:20:200];
649 y0 = [2; 0]; y0d = [0; -2];
650 ng = 1;
651
652 // Simple simulation
653 t = 0:0.003:300;
654 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
655 clf(); plot(yy(1, :), yy(2, :))
656 // Find first point where yy(1) = 0
657 %DAEOPTIONS = list([] , 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
658 [yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
659 plot(yy(1, 1), yy(2, 1), 'r+')
660 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
661
662 // Hot restart for next point
663 t01 = nn(1);
664 [pp, qq] = size(yy);
665 y01 = yy(2:3, qq); y0d1 = yy(3:4, qq);
666 [yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
667 plot(yy(1, 1), yy(2, 1), 'r+')
668 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
669 cd(previous_dir);
670  ]]></programlisting>
671         <scilab:image><![CDATA[
672 code = ['#include <math.h>'
673       'void res22(double *t, double *y, double *yd, double *res, int *ires, double *rpar, int *ipar)'
674       '{res[0] = yd[0] - y[1];'
675       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
676       ' '
677       'void jac22(double *t, double *y, double *yd, double *pd, double *cj, double *rpar, int *ipar)'
678       '{pd[0] = *cj - 0.0;'
679       ' pd[1] =     - (-200.0*y[0]*y[1] - 1.0);'
680       ' pd[2] =     - 1.0;'
681       ' pd[3] = *cj - (100.0*(1.0 - y[0]*y[0]));}'
682       ' '
683       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
684       '{ groot[0] = y[0];}']
685 previous_dir = pwd();
686 cd TMPDIR;
687 mputl(code, 't22.c')
688 ilib_for_link(['res22' 'jac22' 'gr22'], 't22.c', [], 'c', [], 't22loader.sce');
689 exec('t22loader.sce')
690 rtol = [1.d-6; 1.d-6];
691 atol = [1.d-6; 1.d-4];
692 t0 = 0; t = [20:20:200];
693 y0 = [2; 0]; y0d = [0; -2];
694 ng = 1;
695 t = 0:0.003:300;
696 yy = dae([y0, y0d], t0, t, atol, rtol, 'res22', 'jac22');
697 clf(); plot(yy(1, :), yy(2, :))
698 %DAEOPTIONS = list([], 0, [], [], [], 0, [], 1, [], 0, 1, [], [], 1);
699 [yy, nn, hotd] = dae("root2", [y0, y0d], t0, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1');
700 plot(yy(1, 1), yy(2, 1), 'r+')
701 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
702 t01 = nn(1);
703 [pp, qq] = size(yy);
704 y01 = yy(2:3, qq);
705 y0d1 = yy(3:4, qq);
706 [yy, nn, hotd] = dae("root2", [y01, y0d1], t01, 300, atol, rtol, 'res22', 'jac22', ng, 'gr22', 'psol1', 'pjac1', hotd);
707 plot(yy(1, 1), yy(2, 1), 'r+')
708 xstring(yy(1, 1)+0.1, yy(2, 1), string(nn(1)));
709 cd(previous_dir);
710  ]]></scilab:image>
711     </refsection>
712     <refsection role="see also">
713         <title>See Also</title>
714         <simplelist type="inline">
715             <member>
716                 <link linkend="ode">ode</link>
717             </member>
718             <member>
719                 <link linkend="daeoptions">daeoptions</link>
720             </member>
721             <member>
722                 <link linkend="dassl">dassl</link>
723             </member>
724             <member>
725                 <link linkend="dasrt">dasrt</link>
726             </member>
727             <member>
728                 <link linkend="daskr">daskr</link>
729             </member>
730             <member>
731                 <link linkend="impl">impl</link>
732             </member>
733             <member>
734                 <link linkend="fort">fort</link>
735             </member>
736             <member>
737                 <link linkend="link">link</link>
738             </member>
739             <member>
740                 <link linkend="external">external</link>
741             </member>
742         </simplelist>
743     </refsection>
744 </refentry>