c6107af6b8cfe9f194152968ea944de896a5c503
[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) 2008 - INRIA
5  * ...
6  * 
7  * This file must be used under the terms of the CeCILL.
8  * This source file is licensed as described in the file COPYING, which
9  * you should have received as part of this distribution.  The terms
10  * are also available at    
11  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
12  *
13  -->
14 <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">
15     <refnamediv>
16         <refname>dae</refname>
17         <refpurpose>Differential algebraic equations solver</refpurpose>
18     </refnamediv>
19     <refsynopsisdiv>
20         <title>Calling Sequence</title>
21         <synopsis> y=dae(initial,t0,t,res) 
22             [y [,hd]]=dae(initial,t0,t [,rtol, [atol]],res [,jac] [,hd])
23             [y,rd]=dae("root",initial,t0,t,res,ng,surface)
24             [y ,rd [,hd]]=dae("root",initial,t0,t [,rtol, [atol]],res [,jac], ng, surface [,hd])
25         </synopsis>
26     </refsynopsisdiv>
27     <refsection>
28         <title>Arguments</title>
29         <variablelist>
30             <varlistentry>
31                 <term>initial</term>
32                 <listitem>
33                     <para>
34                         a column vector. It may be equal to <literal>x0</literal> or
35                         <literal>[x0;xdot0]</literal>. Where <literal>x0</literal> is the
36                         state value at initial time <literal>t0</literal> and
37                         <literal>xdot0</literal> is the initial state derivative value or an
38                         estimation of it (see below).
39                     </para>
40                 </listitem>
41             </varlistentry>
42             <varlistentry>
43                 <term>t0</term>
44                 <listitem>
45                     <para>a real number, the initial time.</para>
46                 </listitem>
47             </varlistentry>
48             <varlistentry>
49                 <term>t</term>
50                 <listitem>
51                     <para>a real scalar or vector. Gives instants for which you want the
52                         solution. Note that you can get solution at each dae's step point by
53                         setting 
54                         <literal>
55                             <link linkend="daeoptions">%DAEOPTIONS</link>(2)=1
56                         </literal>
57                         .
58                     </para>
59                 </listitem>
60             </varlistentry>
61             <varlistentry>
62                 <term>rtol</term>
63                 <listitem>
64                     <para>
65                         a real scalar or a column vector of same size as
66                         <literal>x0</literal>, the relative error
67                         tolerance of solution. If <literal>rtol</literal> is a vector the
68                         tolerances are specified for each component of the state.
69                     </para>
70                 </listitem>
71             </varlistentry>
72             <varlistentry>
73                 <term>atol</term>
74                 <listitem>
75                     <para>a real scalar or a column vector of same size as
76                         <literal>x0</literal>, the absolute error
77                         tolerance of solution. If <literal>atol</literal> is a vector the
78                         tolerances are specified for each component of the state.
79                     </para>
80                 </listitem>
81             </varlistentry>
82             <varlistentry>
83                 <term>res</term>
84                 <listitem>
85                     <para>
86                         an <link linkend="external" role="" version="">external</link> function computes the value of
87                         <literal>g(t,y,ydot)</literal>. It may be
88                     </para>
89                     <variablelist>
90                         <varlistentry>
91                             <term>a Scilab function</term>
92                             <listitem>
93                                 <para>In this case, its calling sequence must be
94                                     <literal>[r,ires]=res(t,x,xdot)</literal> and
95                                     <literal>res</literal> must return the residue
96                                     <literal>r=g(t,x,xdot)</literal> and error flag
97                                     <literal>ires</literal>.
98                                 </para>
99                                 <para>
100                                     <literal>ires = 0</literal> if
101                                     <literal>res</literal> succeeds to compute
102                                     <literal>r</literal>.
103                                 </para>
104                                 <para>
105                                     <literal>ires = -1</literal> if residue is
106                                     locally not defined for <literal>g(t,x,xdot)</literal>.
107                                 </para>
108                                 <para>
109                                     <literal>ires =-2</literal> if parameters are out of admissible
110                                     range.
111                                 </para>
112                             </listitem>
113                         </varlistentry>
114                         <varlistentry>
115                             <term>a list</term>
116                             <listitem>
117                                 <para>This form of external is used to pass parameters to the
118                                     function. It must be as follows:
119                                 </para>
120                                 <programlisting role="no-scilab-exec"><![CDATA[ 
121 list(res,p1,p2,...)
122  ]]></programlisting>
123                                 <para>where the calling sequence of the function
124                                     <literal>res</literal> is now
125                                 </para>
126                                 <programlisting role="no-scilab-exec"><![CDATA[ 
127 r=res(t,y,ydot,p1,p2,...)
128  ]]></programlisting>
129                                 <para>
130                                     <literal>res</literal> still returns the residual value
131                                     as a function of <literal>(t,x,xdot,x1,x2,...)</literal>, and
132                                     <literal>p1, p2,...</literal> are function parameters.
133                                 </para>
134                             </listitem>
135                         </varlistentry>
136                         <varlistentry>
137                             <term>a character string</term>
138                             <listitem>
139                                 <para>
140                                     it must refer to the name of a C or fortran routine,
141                                     assuming that &lt;<literal>r_name</literal>&gt; is the given name.
142                                 </para>
143                                 <itemizedlist>
144                                     <listitem>
145                                         <para>
146                                             The Fortran calling sequence must be
147                                         </para>
148                                         <para>
149                                             <literal>&lt;r_name&gt;(t,x,xdot,res,ires,rpar,ipar)</literal>
150                                         </para>
151                                         <para>
152                                             <literal>double precision
153                                                 t,x(*),xdot(*),res(*),rpar(*)
154                                             </literal>
155                                         </para>
156                                         <para>
157                                             <literal>integer ires,ipar(*)</literal>
158                                         </para>
159                                     </listitem>
160                                     <listitem>
161                                         <para>The C calling sequence must be</para>
162                                         <para>
163                                             <literal>C2F(&lt;r_name&gt;)(double *t, double *x, double
164                                                 *xdot, double *res, integer *ires, double *rpar, integer
165                                                 *ipar)
166                                             </literal>
167                                         </para>
168                                     </listitem>
169                                 </itemizedlist>
170                                 <para>where</para>
171                                 <itemizedlist>
172                                     <listitem>
173                                         <para>
174                                             <literal>t</literal> is the current time
175                                             value
176                                         </para>
177                                     </listitem>
178                                     <listitem>
179                                         <para>
180                                             <literal>x</literal> the state array
181                                         </para>
182                                     </listitem>
183                                     <listitem>
184                                         <para>
185                                             <literal>xdot</literal> the array of state
186                                             derivatives
187                                         </para>
188                                     </listitem>
189                                     <listitem>
190                                         <para>
191                                             <literal>res</literal> the array of residuals
192                                         </para>
193                                     </listitem>
194                                     <listitem>
195                                         <para>
196                                             <literal>ires</literal> the execution
197                                             indicator
198                                         </para>
199                                     </listitem>
200                                     <listitem>
201                                         <para>
202                                             <literal>rpar</literal> is the array of floating
203                                             point parameter values, needed but cannot be set by the
204                                             <literal>dae</literal> function
205                                         </para>
206                                     </listitem>
207                                     <listitem>
208                                         <para>
209                                             <literal>ipar</literal> is the array of floating
210                                             integer parameter values, needed but cannot be set by the
211                                             <literal>dae</literal> function
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                         an <link linkend="external">external</link> computes the value of
225                         <literal>dg/dx+cj*dg/dxdot</literal> for a given value of parameter
226                         <literal>cj</literal>. It may be
227                     </para>
228                     <variablelist>
229                         <varlistentry>
230                             <term>a Scilab function</term>
231                             <listitem>
232                                 <para>Its calling sequence must be
233                                     <literal>r=jac(t,x,xdot,cj)</literal> and the
234                                     <literal>jac</literal> function must return
235                                     <literal>r=dg(t,x,xdot)/dy+cj*dg(t,x,xdot)/dxdot</literal>
236                                     where <literal>cj</literal> is a real scalar.
237                                 </para>
238                             </listitem>
239                         </varlistentry>
240                         <varlistentry>
241                             <term>a list</term>
242                             <listitem>
243                                 <para>This form of external is used to pass parameters to the
244                                     function. It must be as follows:
245                                 </para>
246                                 <programlisting role="no-scilab-exec"><![CDATA[  
247 list(jac,p1,p2,...)               
248  ]]></programlisting>
249                                 <para>where the calling sequence of the function
250                                     <literal>jac</literal> is now
251                                 </para>
252                                 <programlisting role="no-scilab-exec"><![CDATA[ 
253 r=jac(t,x,xdot,p1,p2,...)      
254  ]]></programlisting>
255                                 <para>
256                                     <literal>jac</literal> still returns
257                                     <literal>dg/dx+cj*dg/dxdot</literal> as a function of
258                                     <literal>(t,x,xdot,cj,p1,p2,...)</literal>.
259                                 </para>
260                             </listitem>
261                         </varlistentry>
262                         <varlistentry>
263                             <term>a character string</term>
264                             <listitem>
265                                 <para>it must refer to the name of a C or fortran routine
266                                     assuming that &lt;j_name&gt; is the given name.
267                                 </para>
268                                 <itemizedlist>
269                                     <listitem>
270                                         <para>
271                                             The Fortran calling sequence must be
272                                         </para>
273                                         <para>
274                                             <literal>&lt;j_name&gt;(t, x, xdot, r, cj, ires,
275                                                 rpar, ipar)
276                                             </literal>
277                                         </para>
278                                         <para>
279                                             double precision <literal>t, x(*), xdot(*), r(*),
280                                                 ci, rpar(*)
281                                             </literal>
282                                         </para>
283                                         <para>
284                                             integer <literal>ires, ipar(*)</literal>
285                                         </para>
286                                     </listitem>
287                                     <listitem>
288                                         <para>The C calling sequence must be</para>
289                                         <para>
290                                             <literal>C2F(&lt;j_name&gt;)(double *t, double *x, double
291                                                 *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)
292                                             </literal>
293                                         </para>
294                                     </listitem>
295                                 </itemizedlist>
296                                 <para>
297                                     where <literal>t, x, xdot, ires, rpar, ipar</literal>
298                                     have similar definition as above, <literal>r</literal> is the results
299                                     array
300                                 </para>
301                             </listitem>
302                         </varlistentry>
303                     </variablelist>
304                 </listitem>
305             </varlistentry>
306             <varlistentry>
307                 <term>surface</term>
308                 <listitem>
309                     <para>
310                         an <link linkend="external">external</link> computes the value of the column vector
311                         <literal>surface(t,x)</literal> with <literal>ng</literal>
312                         components. Each component defines a surface.
313                     </para>
314                     <variablelist>
315                         <varlistentry>
316                             <term>a Scilab function</term>
317                             <listitem>
318                                 <para>Its calling sequence must be
319                                     <literal>r=surface(t,x)</literal>, this function must return a
320                                     vector with <literal>ng</literal> elements.
321                                 </para>
322                             </listitem>
323                         </varlistentry>
324                         <varlistentry>
325                             <term>a list</term>
326                             <listitem>
327                                 <para>
328                                     This form of <link linkend="external">external</link> is used to pass parameters to the
329                                     function. It must be as follows:
330                                 </para>
331                                 <programlisting role="no-scilab-exec"><![CDATA[  
332 list(surface,p1,p2,...)
333  ]]></programlisting>
334                                 <para>where the calling sequence of the function
335                                     <literal>surface</literal> is now
336                                 </para>
337                                 <programlisting role="no-scilab-exec"><![CDATA[ 
338 r=surface(t,x,p1,p2,...)
339  ]]></programlisting>
340                             </listitem>
341                         </varlistentry>
342                         <varlistentry>
343                             <term>a character string</term>
344                             <listitem>
345                                 <para>it must refer to the name of a C or fortran routine.
346                                     Assuming that &lt;s_name&gt; is the given name,
347                                 </para>
348                                 <itemizedlist>
349                                     <listitem>
350                                         <para>
351                                             <literal>The Fortran calling sequence must
352                                                 be
353                                             </literal>
354                                         </para>
355                                         <para>
356                                             <literal>&lt;s_name&gt;(nx, t, x, ng, r, rpar,
357                                                 ipar)
358                                             </literal>
359                                         </para>
360                                         <para>
361                                             <literal>double precision t, x(*), r(*),
362                                                 rpar(*)
363                                             </literal>
364                                         </para>
365                                         <para>
366                                             <literal>integer nx, ng,ipar(*)</literal>
367                                         </para>
368                                     </listitem>
369                                     <listitem>
370                                         <para>The C calling sequence must be</para>
371                                         <para>
372                                             <literal>C2F(&lt;s_name&gt;)(double *t, double *x, double
373                                                 *xdot, double *r, double *cj, integer *ires, double *rpar, integer *ipar)
374                                             </literal>
375                                         </para>
376                                     </listitem>
377                                 </itemizedlist>
378                                 <para>
379                                     where <literal>t, x, rpar, ipar</literal> have similar
380                                     definition as above, <literal>ng</literal> is the number of
381                                     surfaces, <literal>nx</literal> the dimension of the state and
382                                     r is the results array.
383                                 </para>
384                             </listitem>
385                         </varlistentry>
386                     </variablelist>
387                 </listitem>
388             </varlistentry>
389             <varlistentry>
390                 <term>rd</term>
391                 <listitem>
392                     <para>
393                         a vector with two entries <literal>[times num]</literal> where 
394                         <literal>times</literal> is the value of the time at which the
395                         surface is crossed, <literal>num</literal> is the number of the
396                         crossed surface
397                     </para>
398                 </listitem>
399             </varlistentry>
400             <varlistentry>
401                 <term>hd</term>
402                 <listitem>
403                     <para>a real vector, as an output it stores the
404                         <literal>dae</literal> context. It can be used as an input argument
405                         to resume integration (hot restart).
406                     </para>
407                 </listitem>
408             </varlistentry>
409             <varlistentry>
410                 <term>y</term>
411                 <listitem>
412                     <para>
413                         a real matrix. If 
414                         <literal>
415                             <link linkend="daeoptions">%DAEOPTIONS</link>(2)=1
416                         </literal>
417                         ,
418                         each column is the vector <literal>[t;x(t);xdot(t)]</literal> where
419                         <literal>t</literal> is time index for which the solution had been
420                         computed. Else <literal>y</literal> is the vector
421                         <literal>[x(t);xdot(t)]</literal>.
422                     </para>
423                 </listitem>
424             </varlistentry>
425         </variablelist>
426     </refsection>
427     <refsection>
428         <title>Description</title>
429         <para>
430             The <literal>dae</literal> function is a gateway built above the
431             <link linkend="dassl">dassl</link> and <link linkend="dasrt">dasrt</link>
432             function designed for implicit differential equations integration.
433         </para>
434         <programlisting role="no-scilab-exec"><![CDATA[ 
435 g(t,x,xdot)=0
436 x(t0)=x0  and   xdot(t0)=xdot0
437  ]]></programlisting>
438         <para>
439             If <literal>xdot0</literal> is not given in the <emphasis>initial</emphasis>      
440             argument, the <literal>dae</literal> function tries to compute it solving
441             <literal>g(t,x0,xdot0)=0</literal>.
442         </para>
443         <para>
444             if <literal>xdot0</literal> is given in the <emphasis>initial</emphasis>
445             argument it may be either a compatible derivative
446             satisfying <literal>g(t,x0,xdot0)=0</literal> or an approximate value. In the latter case
447             <link linkend="daeoptions">%DAEOPTIONS</link>(7) must be set to 1.
448         </para>
449         <para>Detailed examples using Scilab and C coded externals are given in
450             <literal>modules/differential_equations/tests/unit_tests/dassldasrt.tst</literal>
451         </para>
452     </refsection>
453     <refsection>
454         <title>Examples</title>
455         <programlisting role="example"><![CDATA[ 
456 //Example with Scilab  code
457 function [r,ires]=chemres(t,y,yd)
458     r(1) = -0.04*y(1) + 1d4*y(2)*y(3) - yd(1);
459     r(2) =  0.04*y(1) - 1d4*y(2)*y(3) - 3d7*y(2)*y(2) - yd(2);
460     r(3) =       y(1) +     y(2)      + y(3)-1;
461     ires =  0;
462 endfunction
463 function pd=chemjac(x,y,yd,cj)
464     pd=[-0.04-cj , 1d4*y(3)               , 1d4*y(2);
465          0.04    ,-1d4*y(3)-2*3d7*y(2)-cj ,-1d4*y(2);
466          1       , 1                      , 1       ]
467 endfunction
468
469 x0=[1; 0; 0];
470 xd0=[-0.04; 0.04; 0];
471 t=[1.d-5:0.02:.4, 0.41:.1:4, 40, 400, 4000, 40000, 4d5, 4d6, 4d7, 4d8, 4d9, 4d10];
472
473 y=dae([x0,xd0],0,t,chemres);// returns requested observation time points
474
475 %DAEOPTIONS=list([],1,[],[],[],0,0); // ask  dae mesh points to be returned
476 y=dae([x0,xd0],0,4d10,chemres); // without jacobian
477 y=dae([x0,xd0],0,4d10,chemres,chemjac); // with jacobian
478  ]]></programlisting>
479         <programlisting role="example"><![CDATA[
480
481 //example with C code (C compiler needed) --------------------------------------------------
482 //-1- create the C codes in TMPDIR - Vanderpol equation, implicit form
483 code=['#include <math.h>'
484       'void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)'
485       '{res[0] = yd[0] - y[1];'
486       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
487       ' '
488       'void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)'
489       '{pd[0]=*cj - 0.0;'
490       ' pd[1]=    - (-200.0*y[0]*y[1] - 1.0);'
491       ' pd[2]=    - 1.0;'
492       ' pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}'
493       ' '
494       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
495       '{ groot[0] = y[0];}']
496 cd TMPDIR;
497 mputl(code, 't22.c')
498 //-2- compile and load them
499 ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/t22loader.sce');
500 exec('t22loader.sce')
501 //-3- run
502 rtol=[1.d-6;1.d-6];
503 atol=[1.d-6;1.d-4];
504 t0=0;y0=[2;0];
505 y0d=[0;-2];
506 t=[20:20:200];
507 ng=1;
508 //simple simulation
509 t=0:0.003:300;
510 yy=dae([y0,y0d],t0,t,atol,rtol,'res22','jac22');
511 clf();plot(yy(1,:),yy(2,:))
512 //find first point where yy(1)=0
513 [yy,nn,hotd]=dae("root",[y0,y0d],t0,300,atol,rtol,'res22','jac22',ng,'gr22');
514 plot(yy(1,1),yy(2,1),'r+')
515 xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))
516
517 //hot restart for next point
518 t01=nn(1);
519 [pp,qq]=size(yy);
520 y01=yy(2:3,qq);
521 y0d1=yy(3:4,qq);
522 [yy,nn,hotd]=dae("root",[y01,y0d1],t01,300,atol,rtol,'res22','jac22',ng,'gr22',hotd);
523 plot(yy(1,1),yy(2,1),'r+')
524 xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))
525  ]]></programlisting>
526         <scilab:image><![CDATA[
527 code=['#include <math.h>'
528       'void res22(double *t,double *y,double *yd,double *res,int *ires,double *rpar,int *ipar)'
529       '{res[0] = yd[0] - y[1];'
530       ' res[1] = yd[1] - (100.0*(1.0 - y[0]*y[0])*y[1] - y[0]);}'
531       ' '
532       'void jac22(double *t,double *y,double *yd,double *pd,double *cj,double *rpar,int *ipar)'
533       '{pd[0]=*cj - 0.0;'
534       ' pd[1]=    - (-200.0*y[0]*y[1] - 1.0);'
535       ' pd[2]=    - 1.0;'
536       ' pd[3]=*cj - (100.0*(1.0 - y[0]*y[0]));}'
537       ' '
538       'void gr22(int *neq, double *t, double *y, int *ng, double *groot, double *rpar, int *ipar)'
539       '{ groot[0] = y[0];}']
540 cd TMPDIR;
541 mputl(code, 't22.c')
542 ilib_for_link(['res22' 'jac22' 'gr22'],'t22.c',[],'c',TMPDIR+'/Makefile',TMPDIR+'/t22loader.sce');
543 exec('t22loader.sce')
544 rtol=[1.d-6;1.d-6];
545 atol=[1.d-6;1.d-4];
546 t0=0;y0=[2;0];
547 y0d=[0;-2];
548 t=[20:20:200];
549 ng=1;
550 t=0:0.003:300;
551 yy=dae([y0,y0d],t0,t,atol,rtol,'res22','jac22');
552 clf();plot(yy(1,:),yy(2,:))
553 [yy,nn,hotd]=dae("root",[y0,y0d],t0,300,atol,rtol,'res22','jac22',ng,'gr22');
554 plot(yy(1,1),yy(2,1),'r+')
555 xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))
556 t01=nn(1);
557 [pp,qq]=size(yy);
558 y01=yy(2:3,qq);
559 y0d1=yy(3:4,qq);
560 [yy,nn,hotd]=dae("root",[y01,y0d1],t01,300,atol,rtol,'res22','jac22',ng,'gr22',hotd);
561 plot(yy(1,1),yy(2,1),'r+')
562 xstring(yy(1,1)+0.1,yy(2,1),string(nn(1)))
563  ]]></scilab:image>
564     </refsection>
565     <refsection role="see also">
566         <title>See Also</title>
567         <simplelist type="inline">
568             <member>
569                 <link linkend="ode">ode</link>
570             </member>
571             <member>
572                 <link linkend="daeoptions">daeoptions</link>
573             </member>
574             <member>
575                 <link linkend="dassl">dassl</link>
576             </member>
577             <member>
578                 <link linkend="daskr">daskr</link>
579             </member>
580             <member>
581                 <link linkend="impl">impl</link>
582             </member>
583             <member>
584                 <link linkend="fort">fort</link>
585             </member>
586             <member>
587                 <link linkend="link">link</link>
588             </member>
589             <member>
590                 <link linkend="external">external</link>
591             </member>
592         </simplelist>
593     </refsection>
594 </refentry>